Graphics Plus/Graphics Plus.iso/msdos/animutil/pcmovie/pcmovie.exe

Text File | 1994-01-10 | 379.1 KB | 4,794 lines

//________  JOB
//*
//* APP3.GRAPHICS.PCMOVIE -  8Feb91/mrg
//* from CUSGJES.VIKTOR.SOURCE 8Feb91
//* Installs 77 rtnes in APP1.GRAPHICS.PCMOVIE
//*
//COMPILE  EXEC  FORTC,OPTIONS='VECTOR,DECK,NOOBJECT'
//SYSPUNCH  DD  DISP=(NEW,PASS),UNIT=VIO,DSN=&&OBJIN,
//             SPACE=(TRK,(50,50),RLSE),DCB=OBJECT
//SYSIN     DD *
C GRAPHICS LIBRARY FOR RASTER PLOTS
C WRITTEN FOR IBM 3090 VF - VIKTOR K. DECYK, UCLA
C COPYRIGHT 1990, REGENTS OF THE UNIVERSITY OF CALIFORNIA
C UPDATE: JANUARY 3, 1991
********************************************************************/
*                                                                  */
* This subroutine library was created ad UCLA.                     */
*                                                                  */
* The University of California requires the following disclaimer   */
* concerning all distributed programs:                             */
*                                                                  */
* Although this program material has been tested by its            */
* contributor, no warranty, expressed or implied, is made by the   */
* contributor or the University of California as to the accuracy   */
* and functioning of the program and related program material, nor */
* shall the fact of the distribution constitute any such warranty, */
* and no responsibility is assumed by the contributor or the       */
* University of California, in connection therewith.               */
*                                                                  */
********************************************************************/
*
*****************************************************
* GOPEN -- OPENS GRAPHICS LIBRARY
*****************************************************
*
      SUBROUTINE GOPEN
C THIS SUBROUTINE OPENS GRAPHICS LIBRARY
      COMMON /PLOTCM/ IPLOT, NPLOT, IDTYPE
      IPLOT = 0
      CALL STARTG
      RETURN
      END
*
*****************************************************
* gmopen -- initializes compressed raster device
* USE GMOPEN INSTEAD OF GOPEN FOR THE CONVERTING RASTER IMAGES TO MFE
*****************************************************
*
      SUBROUTINE GMOPEN(IGTYPE,PAL,LPAL)
C THIS SUBROUTINE INITIALIZES COMPRESSED RASTER DEVICE
C FOR MFE FORMAT
C IGTYPE = (1,2,3) = (CGA,EGA,VGA) FORMAT
C PAL = 256 COLOR PALETTE IN RGB FORMAT
C LPAL = LENGTH OF PALETTE (LPAL = 0 MEANS USE DEFAULT PALETTE)
C DEFAULT IS VGA FORMAT
      COMMON /GRDEV/ LX, LY, NBIT, IFRMT, INTRL, IXOR
      common /dithpal/ pal64, npal64
      CHARACTER*1 PAL( lpal*3 + 1 )
      character*1 pal64( 768 )
      CHARACTER*1 C
      DIMENSION LXS(7), LYS(7), NBITS(4)
      SAVE LXS,LYS,NBITS,ISTART
   97 FORMAT (18H PROGRAM EXECUTING)
      DATA LXS /512,640,320,640,720,79,1024/
      DATA LYS /342,480,200,350,384,21,781/
      DATA NBITS /1,2,4,8/
      DATA ISTART /0/
      IF (ISTART.NE.0) GO TO 90
      npal64 = 0
      INTRL  = 0
      IXOR   = 1
      NPAL   = 1
      IFRMT  = 4
      IXOR   = 0

      IF ((IGTYPE.LT.1).OR.(IGTYPE.GT.3)) IGTYPE = 3
      ID = IGTYPE
      IF (ID.EQ.1) then
         NBIT = 2
         INTRL = 1
      else IF (ID.EQ.2) then
         NBIT = 1
      else
         NBIT = 8
         NPAL = 0
      end if

      IF (ID.LT.3) ID = ID + 2
      LX = LXS(ID)
      LY = LYS(ID)
      CALL HEADER(IFRMT,LX,LY,NBIT)

      if( lpal .ne. 0 ) then
         do 71 ijes = 1, lpal*3
            pal64( ijes ) = char(  ichar( pal(ijes) ) / 4   )
71       continue
      end if

      npal64 = lpal
C     write( 6,* ) ' npal = ', npal

      IF (NPAL.EQ.0) then
C        write( 6,* ) ' gmopen calling wrpal '
C        write( 6,* ) ' ifrmt = ', ifrmt
         CALL WRPAL( pal64, npal64, IFRMT )
      end if

      ISTART = 1
   90 WRITE (6,97)
      END
*
*****************************************************
* GRASP1 -- DISPLAYS (X-VX) PHASE SPACE
*****************************************************
*
      SUBROUTINE GRASP1 (PART,LABEL,TIME,VMAX,NX,ITWO,NP,NPX,CI,IRC)
C FOR 1D CODE, THIS SUBROUTINE DISPLAYS (X-VX) PHASE SPACE
      COMMON /PLOTCM/ IPLOT, NPLOT, IDTYPE
      CHARACTER*20 LABEL
      CHARACTER*12 LBLV, LBLU
      CHARACTER*44 LBL
      DIMENSION PART(ITWO,NP)
      SAVE LW
   91 FORMAT (1X,A20,16H PHASE SPACE, T=,F7.2)
      DATA LBLV,LBLU /' VX VERSUS X',' UX VERSUS X'/
      DATA LW /1/
      IRC = 0
      IF (NPLOT.LT.1) GO TO 70
C GET GRAPHICS SIZE PARAMETERS
      CALL GRPARM(1,IRX,IRY,MINX,MINY,LENX,LENY,ICH,ICW,ISLB,ISTCX,ISTCY
     1,NTCX,NTCY,IGSTYL)
C FIND SCALES FOR PLOT
      XMIN = 0.
      XMAX = FLOAT(NX)
      YMIN = -VMAX
      YMAX = VMAX
C FIND LOCATION FOR PLOT
      NPL1 = SQRT(FLOAT(NPLOT-1)) + 0.0001
      NPL = NPL1 + 1
      IY = IPLOT/NPL
      IX = IPLOT - IY*NPL
      MNX = (IX*IRX + MINX)/NPL
      MNY = ((NPL1 - IY)*IRY + MINY)/NPL
      LNX = LENX/NPL
      LNY = LENY/NPL
      JSTCX = ISTCX/NPL
      JSTCY = ISTCY/NPL
      JSLB = ISLB/NPL
      JCH = ICH/NPL
      IF (JCH.LT.1) JCH = 1
      JCW = ICW/NPL
      IF (JCW.LT.1) JCW = 1
      IF (IPLOT.GT.0) GO TO 10
      CALL INITGR(IRX,IRY,JCH,JCW)
      GO TO 20
C DRAW GRID
   10 CALL SELFMP(IRX,IRY)
   20 IC = 7
      CALL TICKS (XMAX,XMIN,YMAX,YMIN,MNX,MNY,LNX,LNY,JCH,JCW,JSLB,JSTCX
     1,JSTCY,NTCX,NTCY,IC,LW,IGSTYL)
C WRITE LABELS
      AT1 = FLOAT(JCH + JCH/3)
      AX = FLOAT(MNX - JSLB)
      AY = FLOAT(MNY - JSTCY) - 2.*AT1
      IT1 = 44
      WRITE (LBL,91) LABEL, TIME
      CALL DRSTRG(LBL,AX,AY,IC,LW,JCW,IT1)
      IT1 = 12
      AY = AY - AT1
      IF (CI.EQ.0) CALL DRSTRG(LBLV,AX,AY,IC,LW,JCW,IT1)
      IF (CI.GT.0) CALL DRSTRG(LBLU,AX,AY,IC,LW,JCW,IT1)
C PLOT GRAPH
      CALL MAPWIN(XMIN,XMAX,YMIN,YMAX,MNX,LNX,MNY,LNY)
      IT2 = 1
      IF (NPX.EQ.0) GO TO 40
      IC = 1
      DO 30 J = 1, NPX
      CALL DRPNTS(PART(1,J),PART(2,J),IT2,IC,LW)
   30 CONTINUE
   40 IF (NPX.GE.NP) GO TO 60
      IC = 4
      IT1 = NPX + 1
      DO 50 J = IT1, NP
      CALL DRPNTS(PART(1,J),PART(2,J),IT2,IC,LW)
   50 CONTINUE
   60 IPLOT = IPLOT + 1
      IF (IPLOT.EQ.NPLOT) IPLOT = 0
      IF (IPLOT.GT.0) GO TO 70
      CALL SGRAPH
      CALL READC(IRC)
   70 RETURN
      END
*
*****************************************************
* GRAF2 -- LINE PLOT OF Y VS X, FILLS MAX DISPLAY REGION
*****************************************************
*
      SUBROUTINE GRAF2 (Y,LABELY,X,LABELX,N,CHR,NCR,IRC)
C THIS SUBROUTINE DOES A LINE PLOT OF Y VERSUS X, WHICH WILL FILL THE
C MAXIMUM ARE OF THE DISPLAY REGION
      COMMON /PLOTCM/ IPLOT, NPLOT, IDTYPE
      CHARACTER*20 LABELY, LABELX
      CHARACTER*48 LBL
      DIMENSION X(N), Y(N)
      CHARACTER*(*) CHR
      SAVE LW,EPS
   91 FORMAT (A20,8H VERSUS ,A20)
      DATA LW /1/
C     DATA EPS /8.0E-14/
      DATA EPS /0./
      IRC = 0
      IF (NPLOT.LT.1) GO TO 40
C GET GRAPHICS SIZE PARAMETERS
      CALL GRPARM(1,IRX,IRY,MINX,MINY,LENX,LENY,ICH,ICW,ISLB,ISTCX,ISTCY
     1,NTCX,NTCY,IGSTYL)
C FIND SCALES FOR PLOT
      XMIN = X(1)
      XMAX = XMIN
      YMIN = Y(1)
      YMAX = YMIN
      DO 10 J = 1, N
      IF (X(J).GT.XMAX) XMAX = X(J)
      IF (X(J).LT.XMIN) XMIN = X(J)
      IF (Y(J).GT.YMAX) YMAX = Y(J)
      IF (Y(J).LT.YMIN) YMIN = Y(J)
   10 CONTINUE
      IF ((XMAX - XMIN).LE.EPS) XMAX = XMIN + 1.
      IF ((YMAX - YMIN).LE.EPS) YMAX = YMIN + 1.
C FIND LOCATION FOR PLOT
      NPL1 = SQRT(FLOAT(NPLOT-1)) + 0.0001
      NPL = NPL1 + 1
      IY = IPLOT/NPL
      IX = IPLOT - IY*NPL
      MNX = (IX*IRX + MINX)/NPL
      MNY = ((NPL1 - IY)*IRY + MINY)/NPL
      LNX = LENX/NPL
      LNY = LENY/NPL
      JSTCX = ISTCX/NPL
      JSTCY = ISTCY/NPL
      JSLB = ISLB/NPL
      JCH = ICH/NPL
      IF (JCH.LT.1) JCH = 1
      JCW = ICW/NPL
      IF (JCW.LT.1) JCW = 1
      IF (IPLOT.GT.0) GO TO 20
      CALL INITGR(IRX,IRY,JCH,JCW)
      GO TO 30
C DRAW GRID
   20 CALL SELFMP(IRX,IRY)
   30 IC = 7
      CALL TICKS (XMAX,XMIN,YMAX,YMIN,MNX,MNY,LNX,LNY,JCH,JCW,JSLB,JSTCX
     1,JSTCY,NTCX,NTCY,IC,LW,IGSTYL)
C WRITE LABELS
      AT1 = FLOAT(JCH + JCH/3)
      AX = FLOAT(MNX - JSLB)
      AY = FLOAT(MNY - JSTCY) - 2.*AT1
      IT1 = 48
      WRITE (LBL,91) LABELY, LABELX
      CALL DRSTRG(LBL,AX,AY,IC,LW,JCW,IT1)
      AY = AY - AT1
      IF (NCR.GT.0) CALL DRSTRG(CHR,AX,AY,IC,LW,JCW,NCR)
C PLOT CURVE
      CALL MAPWIN(XMIN,XMAX,YMIN,YMAX,MNX,LNX,MNY,LNY)
      IC = 1
      CALL DRLINS (X,Y,N,IC,LW)
      IPLOT = IPLOT + 1
      IF (IPLOT.EQ.NPLOT) IPLOT = 0
      IF (IPLOT.GT.0) GO TO 40
      CALL SGRAPH
      CALL READC(IRC)
   40 RETURN
      END
*
*****************************************************
* GRAF1 -- PLOT OF Y VS X WITH SPECIFIED DISPLAY REGION
*****************************************************
*
      SUBROUTINE GRAF1(Y,LABELY,YMAX,YMIN,X,LABELX,XMAX,XMIN,N,CHR,NCR,I
     1RC)
C THIS SUBROUTINE DOES A POINT PLOT OF Y VERSUS X, WITH THE MAXIMUM
C AND MINIMUM VALUES OF THE DISPLAY REGION GIVEN BY YMAX, YMIN, AND
C XMAX, XMIN, RESPECTIVELY.
      COMMON /PLOTCM/ IPLOT, NPLOT, IDTYPE
      CHARACTER*20 LABELY, LABELX
      CHARACTER*48 LBL
      DIMENSION X(N), Y(N)
      CHARACTER*(*) CHR
      SAVE LW,EPS
   91 FORMAT (A20,8H VERSUS ,A20)
      DATA LW /1/
C     DATA EPS /8.0E-14/
      DATA EPS /0./
      IRC = 0
      IF (NPLOT.LT.1) GO TO 30
C GET GRAPHICS SIZE PARAMETERS
      CALL GRPARM(1,IRX,IRY,MINX,MINY,LENX,LENY,ICH,ICW,ISLB,ISTCX,ISTCY
     1,NTCX,NTCY,IGSTYL)
      IF ((XMAX - XMIN).LE.EPS) XMAX = XMIN + 1.
      IF ((YMAX - YMIN).LE.EPS) YMAX = YMIN + 1.
C FIND LOCATION FOR PLOT
      NPL1 = SQRT(FLOAT(NPLOT-1)) + 0.0001
      NPL = NPL1 + 1
      IY = IPLOT/NPL
      IX = IPLOT - IY*NPL
      MNX = (IX*IRX + MINX)/NPL
      MNY = ((NPL1 - IY)*IRY + MINY)/NPL
      LNX = LENX/NPL
      LNY = LENY/NPL
      JSTCX = ISTCX/NPL
      JSTCY = ISTCY/NPL
      JSLB = ISLB/NPL
      JCH = ICH/NPL
      IF (JCH.LT.1) JCH = 1
      JCW = ICW/NPL
      IF (JCW.LT.1) JCW = 1
      IF (IPLOT.GT.0) GO TO 10
      CALL INITGR(IRX,IRY,JCH,JCW)
      GO TO 20
C DRAW GRID
   10 CALL SELFMP(IRX,IRY)
   20 IC = 7
      CALL TICKS (XMAX,XMIN,YMAX,YMIN,MNX,MNY,LNX,LNY,JCH,JCW,JSLB,JSTCX
     1,JSTCY,NTCX,NTCY,IC,LW,IGSTYL)
C WRITE LABELS
      AT1 = FLOAT(JCH + JCH/3)
      AX = FLOAT(MNX - JSLB)
      AY = FLOAT(MNY - JSTCY) - 2.*AT1
      IT1 = 48
      WRITE (LBL,91) LABELY, LABELX
      CALL DRSTRG(LBL,AX,AY,IC,LW,JCW,IT1)
      AY = AY - AT1
      IF (NCR.GT.0) CALL DRSTRG(CHR,AX,AY,IC,LW,JCW,NCR)
C DRAW POINTS
      CALL MAPWIN(XMIN,XMAX,YMIN,YMAX,MNX,LNX,MNY,LNY)
      IC = 1
      CALL DRPNTS (X,Y,N,IC,LW)
      IPLOT = IPLOT + 1
      IF (IPLOT.EQ.NPLOT) IPLOT = 0
      IF (IPLOT.GT.0) GO TO 30
      CALL SGRAPH
      CALL READC(IRC)
   30 RETURN
      END
*
*****************************************************
* GRAF3 -- M LINE PLOTS OF Y VS X
*****************************************************
*
      SUBROUTINE GRAF3 (Y,LABELY,X,LABELX,N,M,NV,CHR,NCR,IRC)
C THIS SUBROUTINE DOES M LINE PLOTS OF SUBARRAYS OF Y VERSUS X, EACH
C PLOT WITH N POINTS, ON A SCALE WHICH WILL FILL THE MAXIMUM AREA OF
C THE DISPLAY REGION.  EACH SUBARRAY IS PLOTTED IN A DIFFERENT COLOR,
C WITH BLUE, RED, YELLOW, CYAN, MAGENTA, AND GREEN USED IN ORDER.
      COMMON /PLOTCM/ IPLOT, NPLOT, IDTYPE
      CHARACTER*20 LABELY, LABELX
      CHARACTER*48 LBL
      DIMENSION ICOLOR(8)
      DIMENSION X(N), Y(NV,M)
      CHARACTER*(*) CHR
      SAVE LW,EPS,ICOLOR
   91 FORMAT (A20,8H VERSUS ,A20)
      DATA LW /1/
C     DATA EPS /8.0E-14/
      DATA EPS /0./
C COLORS ARE: BACKGROUND,BLUE,RED,YELLOW,CYAN,MAGENTA,GREEN,FOREGROUND
      DATA ICOLOR /0,1,4,6,3,5,2,7/
      IRC = 0
      IF (NPLOT.LT.1) GO TO 60
C GET GRAPHICS SIZE PARAMETERS
      CALL GRPARM(1,IRX,IRY,MINX,MINY,LENX,LENY,ICH,ICW,ISLB,ISTCX,ISTCY
     1,NTCX,NTCY,IGSTYL)
C FIND SCALES FOR PLOT
      XMIN = X(1)
      XMAX = XMIN
      YMIN = Y(1,1)
      YMAX = YMIN
      DO 20 J = 1, N
      IF (X(J).GT.XMAX) XMAX = X(J)
      IF (X(J).LT.XMIN) XMIN = X(J)
      DO 10 K = 1, M
      IF (Y(J,K).GT.YMAX) YMAX = Y(J,K)
      IF (Y(J,K).LT.YMIN) YMIN = Y(J,K)
   10 CONTINUE
   20 CONTINUE
      IF ((XMAX - XMIN).LE.EPS) XMAX = XMIN + 1.
      IF ((YMAX - YMIN).LE.EPS) YMAX = YMIN + 1.
C FIND LOCATION FOR PLOT
      NPL1 = SQRT(FLOAT(NPLOT-1)) + 0.0001
      NPL = NPL1 + 1
      IY = IPLOT/NPL
      IX = IPLOT - IY*NPL
      MNX = (IX*IRX + MINX)/NPL
      MNY = ((NPL1 - IY)*IRY + MINY)/NPL
      LNX = LENX/NPL
      LNY = LENY/NPL
      JSTCX = ISTCX/NPL
      JSTCY = ISTCY/NPL
      JSLB = ISLB/NPL
      JCH = ICH/NPL
      IF (JCH.LT.1) JCH = 1
      JCW = ICW/NPL
      IF (JCW.LT.1) JCW = 1
      IF (IPLOT.GT.0) GO TO 30
      CALL INITGR(IRX,IRY,JCH,JCW)
      GO TO 40
C DRAW GRID
   30 CALL SELFMP(IRX,IRY)
   40 IC = 7
      CALL TICKS (XMAX,XMIN,YMAX,YMIN,MNX,MNY,LNX,LNY,JCH,JCW,JSLB,JSTCX
     1,JSTCY,NTCX,NTCY,IC,LW,IGSTYL)
C WRITE LABELS
      AT1 = FLOAT(JCH + JCH/3)
      AX = FLOAT(MNX - JSLB)
      AY = FLOAT(MNY - JSTCY) - 2.*AT1
      IT1 = 48
      WRITE (LBL,91) LABELY, LABELX
      CALL DRSTRG(LBL,AX,AY,IC,LW,JCW,IT1)
      AY = AY - AT1
      IF (NCR.GT.0) CALL DRSTRG(CHR,AX,AY,IC,LW,JCW,NCR)
C PLOT CURVES
      CALL MAPWIN(XMIN,XMAX,YMIN,YMAX,MNX,LNX,MNY,LNY)
      DO 50 K = 1, M
      IT1 = (K - 1)/7
      IC = K - 7*IT1
      CALL DRLINS (X,Y(1,K),N,ICOLOR(IC+1),LW)
   50 CONTINUE
      IPLOT = IPLOT + 1
      IF (IPLOT.EQ.NPLOT) IPLOT = 0
      IF (IPLOT.GT.0) GO TO 60
      CALL SGRAPH
      CALL READC(IRC)
   60 RETURN
      END
*
*****************************************************
* DISP -- PLOTS M SUBARRAYS USING A COMMON SCALE
*****************************************************
*
      SUBROUTINE DISP (F,LABEL,XMIN,XMAX,N,M,NV,ISC,CHR,NCR,IRC)
C THIS SUBROUTINE DISPLAYS M SUBARRAYS OF THE ARRAY F, EACH PLOT WITH N
C POINTS, ON A COMMON SCALE GIVEN BY YMAX = 2**ISC, YMIN = -2**ISC,
C VERSUS A LINEAR FUNCTION IN X, WHERE XMIN < X < XMAX.  IF ABS(ISC) >
C 116, THEN THE PROGRAM FINDS THE MINIMUM VALUE OF ISC WHICH WILL CON-
C TAIN THE PLOTS.  EACH SUBARRAY IS PLOTTED IN A DIFFERENT COLOR,
C WITH BLUE, RED, YELLOW, CYAN, MAGENTA, AND GREEN USED IN ORDER.
      COMMON /PLOTCM/ IPLOT, NPLOT, IDTYPE
      CHARACTER*20 LABEL
      CHARACTER*36 LBL
      DIMENSION F(NV,M)
      CHARACTER*(*) CHR
      DIMENSION ICOLOR(8)
      DIMENSION X(2)
      SAVE LW,DV,ICOLOR
   91 FORMAT (A20,8H, SCALE=,I8)
      DATA LW /1/
      DATA DV /2./
C COLORS ARE: BACKGROUND,BLUE,RED,YELLOW,CYAN,MAGENTA,GREEN,FOREGROUND
      DATA ICOLOR /0,1,4,6,3,5,2,7/
      IRC = 0
      IF (NPLOT.LT.1) GO TO 90
C GET GRAPHICS SIZE PARAMETERS
      CALL GRPARM(1,IRX,IRY,MINX,MINY,LENX,LENY,ICH,ICW,ISLB,ISTCX,ISTCY
     1,NTCX,NTCY,IGSTYL)
C FIND SCALES FOR PLOT
      IS = ISC
      IF (IABS(IS).LE.116) GO TO 30
      FMAX = ABS(F(1,1))
      DO 20 J = 1, M
      DO 10 I = 1, N
      AT1 = ABS(F(I,J))
      IF (AT1.GT.FMAX) FMAX = AT1
   10 CONTINUE
   20 CONTINUE
      IF (FMAX.EQ.0.) FMAX = 1.0E-35
      IS = ALOG(FMAX)/ALOG(DV)
      IF (FMAX.GE.1.) IS = IS + 1
   30 YMAX = DV**IS
      YMIN = -YMAX
C FIND LOCATION FOR PLOT
      NPL1 = SQRT(FLOAT(NPLOT-1)) + 0.0001
      NPL = NPL1 + 1
      IY = IPLOT/NPL
      IX = IPLOT - IY*NPL
      MNX = (IX*IRX + MINX)/NPL
      MNY = ((NPL1 - IY)*IRY + MINY)/NPL
      LNX = LENX/NPL
      LNY = LENY/NPL
      JSTCX = ISTCX/NPL
      JSTCY = ISTCY/NPL
      JSLB = ISLB/NPL
      JCH = ICH/NPL
      IF (JCH.LT.1) JCH = 1
      JCW = ICW/NPL
      IF (JCW.LT.1) JCW = 1
      IF (IPLOT.GT.0) GO TO 40
      CALL INITGR(IRX,IRY,JCH,JCW)
      GO TO 50
C DRAW GRID
   40 CALL SELFMP(IRX,IRY)
   50 IC = 7
      CALL TICKS (XMAX,XMIN,YMAX,YMIN,MNX,MNY,LNX,LNY,JCH,JCW,JSLB,JSTCX
     1,JSTCY,NTCX,NTCY,IC,LW,IGSTYL)
C WRITE LABELS
      AT1 = FLOAT(JCH + JCH/3)
      AX = FLOAT(MNX - JSLB)
      AY = FLOAT(MNY - JSTCY) - 2.*AT1
      IT1 = 36
      WRITE (LBL,91) LABEL, IS
      CALL DRSTRG(LBL,AX,AY,IC,LW,JCW,IT1)
      AY = AY - AT1
      IF (NCR.GT.0) CALL DRSTRG(CHR,AX,AY,IC,LW,JCW,NCR)
C DRAW CURVES
      CALL MAPWIN(XMIN,XMAX,YMIN,YMAX,MNX,LNX,MNY,LNY)
      IF (N.LT.2) GO TO 80
      IT2 = 2
      N1 = N - 1
      DX = (XMAX - XMIN)/FLOAT(N1)
      DO 70 K = 1, M
      IT1 = (K - 1)/7
      IC = K - 7*IT1
      IT3 = ICOLOR(IC+1)
      X(2) = XMIN
      DO 60 I = 1, N1
      X(1) = X(2)
      X(2) = XMIN + DX*FLOAT(I)
      CALL DRLINS (X,F(I,K),IT2,IT3,LW)
   60 CONTINUE
   70 CONTINUE
   80 IPLOT = IPLOT + 1
      IF (IPLOT.EQ.NPLOT) IPLOT = 0
      IF (IPLOT.GT.0) GO TO 90
      CALL SGRAPH
      CALL READC(IRC)
   90 RETURN
      END
*
*****************************************************
* DISP1 -- DISPLAYS M SUBARRAYS OF ARRAY F
*****************************************************
*
      SUBROUTINE DISP1 (F,LABEL,XMIN,XMAX,N,M,NV,ISC,IST,CHR,NCR,IRC)
C THIS SUBROUTINE DISPLAYS M SUBARRAYS OF THE ARRAY F, EACH PLOT WITH N
C POINTS, VERSUS A LINEAR FUNCTION IN X, WHERE XMIN < X < XMAX.
C THE PLOTS HAVE A COMMON SCALE IN Y GIVEN BY YMAX AND YMIN.
C IF IST = 0, THEN YMAX = 2**ISC AND YMIN = -2**ISC.
C IF IST > 0, THEN YMAX = 2**ISC AND YMIN = 0.
C IF IST < 0, THEN YMAX = 0 AND YMIN = -2**ISC.
C IF ABS(ISC) > 116, THEN THE PROGRAM FINDS THE MINIMUM VALUE OF ISC
C WHICH WILL CONTAIN THE PLOTS.  EACH SUBARRAY IS PLOTTED IN A DIFFERENT
C COLOR, WITH BLUE, RED, YELLOW, CYAN, MAGENTA, AND GREEN USED IN ORDER.
      COMMON /PLOTCM/ IPLOT, NPLOT, IDTYPE
      CHARACTER*20 LABEL
      CHARACTER*36 LBL
      DIMENSION F(NV,M)
      CHARACTER*(*) CHR
      DIMENSION ICOLOR(8)
      DIMENSION X(2)
      SAVE LW,DV,ICOLOR
   91 FORMAT (A20,8H, SCALE=,I8)
      DATA LW /1/
      DATA DV /2./
C COLORS ARE: BACKGROUND,BLUE,RED,YELLOW,CYAN,MAGENTA,GREEN,FOREGROUND
      DATA ICOLOR /0,1,4,6,3,5,2,7/
      IRC = 0
      IF (NPLOT.LT.1) GO TO 90
C GET GRAPHICS SIZE PARAMETERS
      CALL GRPARM(1,IRX,IRY,MINX,MINY,LENX,LENY,ICH,ICW,ISLB,ISTCX,ISTCY
     1,NTCX,NTCY,IGSTYL)
C FIND SCALES FOR PLOT
      IS = ISC
      IF (IABS(IS).LE.116) GO TO 30
      FMAX = ABS(F(1,1))
      DO 20 J = 1, M
      DO 10 I = 1, N
      AT1 = ABS(F(I,J))
      IF (AT1.GT.FMAX) FMAX = AT1
   10 CONTINUE
   20 CONTINUE
      IF (FMAX.EQ.0.) FMAX = 1.0E-35
      IS = ALOG(FMAX)/ALOG(DV)
      IF (FMAX.GE.1.) IS = IS + 1
   30 YMAX = DV**IS
      YMIN = -YMAX
      IF (IST.GT.0) YMIN = 0.
      IF (IST.LT.0) YMAX = 0.
C FIND LOCATION FOR PLOT
      NPL1 = SQRT(FLOAT(NPLOT-1)) + 0.0001
      NPL = NPL1 + 1
      IY = IPLOT/NPL
      IX = IPLOT - IY*NPL
      MNX = (IX*IRX + MINX)/NPL
      MNY = ((NPL1 - IY)*IRY + MINY)/NPL
      LNX = LENX/NPL
      LNY = LENY/NPL
      JSTCX = ISTCX/NPL
      JSTCY = ISTCY/NPL
      JSLB = ISLB/NPL
      JCH = ICH/NPL
      IF (JCH.LT.1) JCH = 1
      JCW = ICW/NPL
      IF (JCW.LT.1) JCW = 1
      IF (IPLOT.GT.0) GO TO 40
      CALL INITGR(IRX,IRY,JCH,JCW)
      GO TO 50
C DRAW GRID
   40 CALL SELFMP(IRX,IRY)
   50 IC = 7
      CALL TICKS (XMAX,XMIN,YMAX,YMIN,MNX,MNY,LNX,LNY,JCH,JCW,JSLB,JSTCX
     1,JSTCY,NTCX,NTCY,IC,LW,IGSTYL)
C WRITE LABELS
      AT1 = FLOAT(JCH + JCH/3)
      AX = FLOAT(MNX - JSLB)
      AY = FLOAT(MNY - JSTCY) - 2.*AT1
      IT1 = 36
      WRITE (LBL,91) LABEL, IS
      CALL DRSTRG(LBL,AX,AY,IC,LW,JCW,IT1)
      AY = AY - AT1
      IF (NCR.GT.0) CALL DRSTRG(CHR,AX,AY,IC,LW,JCW,NCR)
C DRAW CURVES
      CALL MAPWIN(XMIN,XMAX,YMIN,YMAX,MNX,LNX,MNY,LNY)
      IF (N.LT.2) GO TO 80
      IT2 = 2
      N1 = N - 1
      DX = (XMAX - XMIN)/FLOAT(N1)
      DO 70 K = 1, M
      IT1 = (K - 1)/7
      IC = K - 7*IT1
      IT3 = ICOLOR(IC+1)
      X(2) = XMIN
      DO 60 I = 1, N1
      X(1) = X(2)
      X(2) = XMIN + DX*FLOAT(I)
      CALL DRLINS (X,F(I,K),IT2,IT3,LW)
   60 CONTINUE
   70 CONTINUE
   80 IPLOT = IPLOT + 1
      IF (IPLOT.EQ.NPLOT) IPLOT = 0
      IF (IPLOT.GT.0) GO TO 90
      CALL SGRAPH
      CALL READC(IRC)
   90 RETURN
      END
*
*****************************************************
* DISP2 -- PLOTS TWO SUBARRAYS OF F
*****************************************************
*
      SUBROUTINE DISP2 (F,LABEL,XV,XMIN,XMAX,N,M,NV,ISC,CHR,NCR,IRC)
C THIS SUBROUTINE DISPLAYS TWO SUBARRAYS OF THE ARRAY F, BOTH PLOTS WITH
C N POINTS, ON A COMMON SCALE GIVEN BY YMAX = 2**ISC, YMIN =-2**ISC,
C VERSUS A LINEAR FUNCTION IN X, WHERE XMIN < X < XMAX.
C IF ABS(ISC) > 116, THEN THE PROGRAM FINDS THE MINIMUM VALUE OF ISC
C WHICH WILL CONTAIN THE PLOTS.  THE FIRST SUBARRAY IS DRAWN AS A LINE
C PLOT, AND THE SECOND AS SMALL CLOSED CIRCLES.  IN ADDITION, A VERTICAL
C LINE AT LOCATION X = XV IS DRAWN.  THE FIRST SUBARRAY IS PLOTTED IN
C BLUE AND THE SECOND IN RED.
      COMMON /PLOTCM/ IPLOT, NPLOT, IDTYPE
      CHARACTER*1 CHS
      CHARACTER*20 LABEL
      CHARACTER*36 LBL
      DIMENSION F(NV,M)
      CHARACTER*(*) CHR
      DIMENSION X(2), Y(2)
      SAVE LW,DV,CHS
   91 FORMAT (A20,8H, SCALE=,I8)
      DATA LW /1/
      DATA DV /2./
      DATA CHS /'o'/
      IRC = 0
      IF (NPLOT.LT.1) GO TO 90
C GET GRAPHICS SIZE PARAMETERS
      CALL GRPARM(1,IRX,IRY,MINX,MINY,LENX,LENY,ICH,ICW,ISLB,ISTCX,ISTCY
     1,NTCX,NTCY,IGSTYL)
C FIND SCALES FOR PLOT
      IS = ISC
      IF (IABS(IS).LE.116) GO TO 30
      FMAX = ABS(F(1,1))
      DO 20 J = 1, M
      DO 10 I = 1, N
      AT1 = ABS(F(I,J))
      IF (AT1.GT.FMAX) FMAX = AT1
   10 CONTINUE
   20 CONTINUE
      IF (FMAX.EQ.0.) FMAX = 1.0E-35
      IS = ALOG(FMAX)/ALOG(DV)
      IF (FMAX.GE.1.) IS = IS + 1
   30 YMAX = DV**IS
      YMIN = -YMAX
C FIND LOCATION FOR PLOT
      NPL1 = SQRT(FLOAT(NPLOT-1)) + 0.0001
      NPL = NPL1 + 1
      IY = IPLOT/NPL
      IX = IPLOT - IY*NPL
      MNX = (IX*IRX + MINX)/NPL
      MNY = ((NPL1 - IY)*IRY + MINY)/NPL
      LNX = LENX/NPL
      LNY = LENY/NPL
      JSTCX = ISTCX/NPL
      JSTCY = ISTCY/NPL
      JSLB = ISLB/NPL
      JCH = ICH/NPL
      IF (JCH.LT.1) JCH = 1
      JCW = ICW/NPL
      IF (JCW.LT.1) JCW = 1
      IF (IPLOT.GT.0) GO TO 40
      CALL INITGR(IRX,IRY,JCH,JCW)
      GO TO 50
C DRAW GRID
   40 CALL SELFMP(IRX,IRY)
   50 IC = 7
      CALL TICKS (XMAX,XMIN,YMAX,YMIN,MNX,MNY,LNX,LNY,JCH,JCW,JSLB,JSTCX
     1,JSTCY,NTCX,NTCY,IC,LW,IGSTYL)
C WRITE LABELS
      AT1 = FLOAT(JCH + JCH/3)
      AX = FLOAT(MNX - JSLB)
      AY = FLOAT(MNY - JSTCY) - 2.*AT1
      IT1 = 36
      WRITE (LBL,91) LABEL, IS
      CALL DRSTRG(LBL,AX,AY,IC,LW,JCW,IT1)
      AY = AY - AT1
      IF (NCR.GT.0) CALL DRSTRG(CHR,AX,AY,IC,LW,JCW,NCR)
      CALL MAPWIN(XMIN,XMAX,YMIN,YMAX,MNX,LNX,MNY,LNY)
C SPECIAL CASE OF VERTICAL LINE
      X(1) = XV
      Y(1) = YMIN
      X(2) = XV
      Y(2) = YMAX
      CALL DRLINS (X,Y,2,IC,LW)
C MULTIPLE CURVES
      IF (N.LT.2) GO TO 80
      IT1 = 2
      N1 = N - 1
      DX = (XMAX - XMIN)/FLOAT(N1)
C FIRST CURVE IS SOLID LINE
      IC = 1
      X(2) = XMIN
      DO 60 I = 1, N1
      X(1) = X(2)
      X(2) = XMIN + DX*FLOAT(I)
      CALL DRLINS (X,F(I,1),IT1,IC,LW)
   60 CONTINUE
C SECOND CURVE IS POINTS
      IF (M.EQ.1) GO TO 80
      IC = 4
      IT1 = 1
      DXH = .3*FLOAT(JCW)*(XMAX - XMIN)/FLOAT(LNX)
      DYH = .25*FLOAT(JCH)*(YMAX - YMIN)/FLOAT(LNY)
      DO 70 I = 1, N
      AX = XMIN + DX*FLOAT(I - 1) - DXH
      AY = F(I,2) - DYH
      CALL DRSTRG(CHS,AX,AY,IC,LW,JCW,IT1)
   70 CONTINUE
   80 IPLOT = IPLOT + 1
      IF (IPLOT.EQ.NPLOT) IPLOT = 0
      IF (IPLOT.GT.0) GO TO 90
      CALL SGRAPH
      CALL READC(IRC)
   90 RETURN
      END
*
*****************************************************
* CONTUR -- CONTOUR PLOT OF FUNCTION F
*****************************************************
*
      SUBROUTINE CONTUR (F,LINK,LABEL,NX,NY,NC,NXV,CHR,NCR,IRC)
C CONTUR DOES A CONTOUR PLOT OF THE FUNCTION F, FOR VALUES OF THE FIRST
C INDEX IN THE RANGE 1 < I < N, AND OF THE SECOND INDEX 1 < J < M.  NC
C CONTOUR INTERVALS ARE CHOSEN, SPACED EQUALLY BETWEEN THE MAXIMUM AND
C MINIMUM VALUES OF F.  SEVEN COLORS ARE USED TO PLOT THE CONTOURS.
C RED, MAGENTA, YELLOW, FOREGROUND, CYAN, GREEN, AND BLUE ARE USED IN
C GOING FROM HIGHEST TO LOWEST VALUES.
C ISTYLE = (0,1) = CONTOUR PLOT (FILLS AREA,PRESERVES ASPECT RATIO)
      COMMON /PLOTCM/ IPLOT, NPLOT, IDTYPE
      LOGICAL*1 LINK(2,NX,NY)
      CHARACTER*20 LABEL
      CHARACTER*57 LBL
      DIMENSION F(NXV,NY)
      CHARACTER*(*) CHR
      DIMENSION X(2), Y(2), C(2)
   91 FORMAT (A20,5H MAX=,E10.3,5H MIN=,E10.3,4H NC=,I3)
      SAVE LW,ISTYLE,ZERO,ONE
      DATA LW,ISTYLE /1,1/
      DATA ZERO,ONE /0.,1./
      IRC = 0
      IF (NPLOT.LT.1) GO TO 70
      ANX = FLOAT(NX)
      ANY = FLOAT(NY)
C GET GRAPHICS SIZE PARAMETERS
      CALL GRPARM(2,IRX,IRY,MINX,MINY,LENX,LENY,ICH,ICW,ISLB,ISTCX,ISTCY
     1,NTCX,NTCY,IGSTYL)
C FIND SCALES FOR PLOT
      FMIN = F(1,1)
      FMAX = FMIN
      DO 20 K = 1, NY
      DO 10 J = 1, NX
      IF (F(J,K).GT.FMAX) FMAX = F(J,K)
      IF (F(J,K).LT.FMIN) FMIN = F(J,K)
   10 CONTINUE
   20 CONTINUE
C FIND LOCATION FOR PLOT
      NPL1 = SQRT(FLOAT(NPLOT-1)) + 0.0001
      NPL = NPL1 + 1
      IY = IPLOT/NPL
      IX = IPLOT - IY*NPL
      MNX = (IX*IRX + MINX)/NPL
      MNY = ((NPL1 - IY)*IRY + MINY)/NPL
      LNX = LENX/NPL
      LNY = LENY/NPL
      JSLB = ISLB/NPL
      JCH = ICH/NPL
      IF (JCH.LT.1) JCH = 1
      JCW = ICW/NPL
      IF (JCW.LT.1) JCW = 1
      IF (IPLOT.GT.0) GO TO 30
      CALL INITGR(IRX,IRY,JCH,JCW)
      GO TO 40
C DRAW GRID
   30 CALL SELFMP(IRX,IRY)
   40 IC = 7
      MNX0 = MNX
      MNY0 = MNY
      IF (ISTYLE.EQ.0) GO TO 50
      MXX = MNX + LNX
      MXY = MNY + LNY
      IF (NY.GT.NX) LNX = FLOAT(LNX)*(ANX/ANY) + .5
      IF (NX.GT.NY) LNY = FLOAT(LNY)*(ANY/ANX) + .5
      MNX = MXX - LNX
      MNY = MXY - LNY
C DRAW BOX
   50 CALL BOX (MNX,MNY,LNX,LNY,JCH,JCW,IC,LW,IGSTYL)
C WRITE LABELS
      AT1 = FLOAT(JCH + JCH/3)
      AX = FLOAT(MNX0 - JSLB)
      AY = FLOAT(MNY0) - AT1
      IT1 = 57
      WRITE (LBL,91) LABEL, FMAX, FMIN, NC
      CALL DRSTRG(LBL,AX,AY,IC,LW,JCW,IT1)
      AY = AY - AT1
      IF (NCR.GT.0) CALL DRSTRG(CHR,AX,AY,IC,LW,JCW,NCR)
C DRAW CONTOURS
      IF (FMAX.EQ.FMIN) GO TO 60
      CALL MAPWIN(ZERO,ANX,ZERO,ANY,MNX,LNX,MNY,LNY)
C FIX STARTING VALUE AND INCREMENT FOR SPATIAL COORDINATES
      X(1) = ZERO
      Y(1) = ZERO
      X(2) = ONE
      Y(2) = ONE
C FIX FIRST CONTOUR LEVEL AND CONTOUR INTERVAL
      C(2) = (FMAX - FMIN)/FLOAT(NC)
      C(1) = FMIN + .5*C(2)
      CALL CONTRU (X,Y,F,C,NX,NY,NC,LINK,NXV,LW)
   60 IPLOT = IPLOT + 1
      IF (IPLOT.EQ.NPLOT) IPLOT = 0
      IF (IPLOT.GT.0) GO TO 70
      CALL SGRAPH
      CALL READC(IRC)
   70 RETURN
      END
*
*****************************************************
* DISPCN -- CONTOUR PLOT OF FUNCTION F
*****************************************************
*
      SUBROUTINE DISPCN (F,LINK,LABEL,NX,NY,NC,NXV,ISC,CHR,NCR,IRC)
C DISPCN DOES A CONTOUR PLOT OF THE FUNCTION F, FOR VALUES OF THE FIRST
C INDEX IN THE RANGE 1 < I < N, AND OF THE SECOND INDEX 1 < J < M.  NC
C CONTOUR INTERVALS ARE CHOSEN, SPACED EQUALLY BETWEEN FMAX = 2**ISC AND
C FMIN = -2**ISC.  IF ABS(ISC) > 116, THEN THE PROGRAM FINDS THE MINIMUM
C VALUE OF ISC WHICH WILL CONTAIN THE FUNCTION VALUES.  SEVEN COLORS ARE
C USED TO PLOT THE CONTOURS.  RED, MAGENTA, YELLOW, FOREGROUND, CYAN,
C GREEN, AND BLUE ARE USED IN GOING FROM HIGHEST TO LOWEST VALUES.
C ISTYLE = (0,1) = CONTOUR PLOT (FILLS AREA,PRESERVES ASPECT RATIO)
      COMMON /PLOTCM/ IPLOT, NPLOT, IDTYPE
      LOGICAL*1 LINK(2,NX,NY)
      CHARACTER*20 LABEL
      CHARACTER*40 LBL
      DIMENSION F(NXV,NY)
      CHARACTER*(*) CHR
      DIMENSION X(2), Y(2), C(2)
   91 FORMAT (A20,8H, SCALE=,I5,4H NC=,I3)
      SAVE LW,ISTYLE,DV,ZERO,ONE
      DATA LW,ISTYLE /1,1/
      DATA DV /2./
      DATA ZERO,ONE /0.,1./
      IRC = 0
      IF (NPLOT.LT.1) GO TO 70
      ANX = FLOAT(NX)
      ANY = FLOAT(NY)
C GET GRAPHICS SIZE PARAMETERS
      CALL GRPARM(2,IRX,IRY,MINX,MINY,LENX,LENY,ICH,ICW,ISLB,ISTCX,ISTCY
     1,NTCX,NTCY,IGSTYL)
C FIND SCALES FOR PLOT
      IS = ISC
      IF (IABS(IS).LE.116) GO TO 30
      FMAX = ABS(F(1,1))
      DO 20 K = 1, NY
      DO 10 J = 1, NX
      AT1 = ABS(F(J,K))
      IF (AT1.GT.FMAX) FMAX = AT1
   10 CONTINUE
   20 CONTINUE
      IF (FMAX.EQ.0.) FMAX = 1.0E-35
      IS = ALOG(FMAX)/ALOG(DV)
      IF (FMAX.GE.1.) IS = IS + 1
   30 FMAX = DV**IS
      FMIN = -FMAX
C FIND LOCATION FOR PLOT
      NPL1 = SQRT(FLOAT(NPLOT-1)) + 0.0001
      NPL = NPL1 + 1
      IY = IPLOT/NPL
      IX = IPLOT - IY*NPL
      MNX = (IX*IRX + MINX)/NPL
      MNY = ((NPL1 - IY)*IRY + MINY)/NPL
      LNX = LENX/NPL
      LNY = LENY/NPL
      JSLB = ISLB/NPL
      JCH = ICH/NPL
      IF (JCH.LT.1) JCH = 1
      JCW = ICW/NPL
      IF (JCW.LT.1) JCW = 1
      IF (IPLOT.GT.0) GO TO 40
      CALL INITGR(IRX,IRY,JCH,JCW)
      GO TO 50
C DRAW GRID
   40 CALL SELFMP(IRX,IRY)
   50 IC = 7
      MNX0 = MNX
      MNY0 = MNY
      IF (ISTYLE.EQ.0) GO TO 60
      MXX = MNX + LNX
      MXY = MNY + LNY
      IF (NY.GT.NX) LNX = FLOAT(LNX)*(ANX/ANY) + .5
      IF (NX.GT.NY) LNY = FLOAT(LNY)*(ANY/ANX) + .5
      MNX = MXX - LNX
      MNY = MXY - LNY
C DRAW BOX
   60 CALL BOX (MNX,MNY,LNX,LNY,JCH,JCW,IC,LW,IGSTYL)
C WRITE LABELS
      AT1 = FLOAT(JCH + JCH/3)
      AX = FLOAT(MNX0 - JSLB)
      AY = FLOAT(MNY0) - AT1
      IT1 = 40
      WRITE (LBL,91) LABEL, IS, NC
      CALL DRSTRG(LBL,AX,AY,IC,LW,JCW,IT1)
      AY = AY - AT1
      IF (NCR.GT.0) CALL DRSTRG(CHR,AX,AY,IC,LW,JCW,NCR)
C DRAW CONTOURS
      CALL MAPWIN(ZERO,ANX,ZERO,ANY,MNX,LNX,MNY,LNY)
C FIX STARTING VALUE AND INCREMENT FOR SPATIAL COORDINATES
      X(1) = ZERO
      Y(1) = ZERO
      X(2) = ONE
      Y(2) = ONE
C FIX FIRST CONTOUR LEVEL AND CONTOUR INTERVAL
      C(2) = (FMAX - FMIN)/FLOAT(NC)
      C(1) = FMIN + .5*C(2)
      CALL CONTRU (X,Y,F,C,NX,NY,NC,LINK,NXV,LW)
      IPLOT = IPLOT + 1
      IF (IPLOT.EQ.NPLOT) IPLOT = 0
      IF (IPLOT.GT.0) GO TO 70
      CALL SGRAPH
      CALL READC(IRC)
   70 RETURN
      END
*
*****************************************************
* DSPCN2 -- CONTOUR PLOT OF FUNCTION F
*****************************************************
*
      SUBROUTINE DSPCN2 (F,LINK,LABEL,XV,YV,NX,NY,NC,NXV,ISC,CHR,NCR,IRC
     1)
C DSPCN2 DOES A CONTOUR PLOT OF THE FUNCTION F, FOR VALUES OF THE FIRST
C INDEX IN THE RANGE 1 < I < N, AND OF THE SECOND INDEX 1 < J < M.  NC
C CONTOUR INTERVALS ARE CHOSEN, SPACED EQUALLY BETWEEN FMAX = 2**ISC AND
C FMIN = -2**ISC.  IF ABS(ISC) > 116, THEN THE PROGRAM FINDS THE MINIMUM
C VALUE OF ISC WHICH WILL CONTAIN THE FUNCTION VALUES.  IN ADDITION, A
C PAIR OF VERTICAL AND HORIZONTAL LINES ARE DRAWN AT X = XV AND Y = YV,
C WHERE THE UNITS ARE 0 < X < N AND 0 < Y < M.  SEVEN COLORS ARE USED TO
C PLOT THE CONTOURS.  RED, MAGENTA, YELLOW, FOREGROUND, CYAN, GREEN, AND
C BLUE ARE USED IN GOING FROM HIGHEST TO LOWEST VALUES.
C ISTYLE = (0,1) = CONTOUR PLOT (FILLS AREA,PRESERVES ASPECT RATIO)
      COMMON /PLOTCM/ IPLOT, NPLOT, IDTYPE
      LOGICAL*1 LINK(2,NX,NY)
      CHARACTER*20 LABEL
      CHARACTER*40 LBL
      DIMENSION F(NXV,NY)
      CHARACTER*(*) CHR
      DIMENSION X(2), Y(2), C(2)
   91 FORMAT (A20,8H, SCALE=,I5,4H NC=,I3)
      SAVE LW,ISTYLE,DV,ZERO,ONE
      DATA LW,ISTYLE /1,1/
      DATA DV /2./
      DATA ZERO,ONE /0.,1./
      IRC = 0
      IF (NPLOT.LT.1) GO TO 70
      ANX = FLOAT(NX)
      ANY = FLOAT(NY)
C GET GRAPHICS SIZE PARAMETERS
      CALL GRPARM(2,IRX,IRY,MINX,MINY,LENX,LENY,ICH,ICW,ISLB,ISTCX,ISTCY
     1,NTCX,NTCY,IGSTYL)
C FIND SCALES FOR PLOT
      IS = ISC
      IF (IABS(IS).LE.116) GO TO 30
      FMAX = ABS(F(1,1))
      DO 20 K = 1, NY
      DO 10 J = 1, NX
      AT1 = ABS(F(J,K))
      IF (AT1.GT.FMAX) FMAX = AT1
   10 CONTINUE
   20 CONTINUE
      IF (FMAX.EQ.0.) FMAX = 1.0E-35
      IS = ALOG(FMAX)/ALOG(DV)
      IF (FMAX.GE.1.) IS = IS + 1
   30 FMAX = DV**IS
      FMIN = -FMAX
C FIND LOCATION FOR PLOT
      NPL1 = SQRT(FLOAT(NPLOT-1)) + 0.0001
      NPL = NPL1 + 1
      IY = IPLOT/NPL
      IX = IPLOT - IY*NPL
      MNX = (IX*IRX + MINX)/NPL
      MNY = ((NPL1 - IY)*IRY + MINY)/NPL
      LNX = LENX/NPL
      LNY = LENY/NPL
      JSLB = ISLB/NPL
      JCH = ICH/NPL
      IF (JCH.LT.1) JCH = 1
      JCW = ICW/NPL
      IF (JCW.LT.1) JCW = 1
      IF (IPLOT.GT.0) GO TO 40
      CALL INITGR(IRX,IRY,JCH,JCW)
      GO TO 50
C DRAW GRID
   40 CALL SELFMP(IRX,IRY)
   50 IC = 7
      MNX0 = MNX
      MNY0 = MNY
      IF (ISTYLE.EQ.0) GO TO 60
      MXX = MNX + LNX
      MXY = MNY + LNY
      IF (NY.GT.NX) LNX = FLOAT(LNX)*(ANX/ANY) + .5
      IF (NX.GT.NY) LNY = FLOAT(LNY)*(ANY/ANX) + .5
      MNX = MXX - LNX
      MNY = MXY - LNY
C DRAW BOX
   60 CALL BOX (MNX,MNY,LNX,LNY,JCH,JCW,IC,LW,IGSTYL)
C WRITE LABELS
      AT1 = FLOAT(JCH + JCH/3)
      AX = FLOAT(MNX0 - JSLB)
      AY = FLOAT(MNY0) - AT1
      IT1 = 40
      WRITE (LBL,91) LABEL, IS, NC
      CALL DRSTRG(LBL,AX,AY,IC,LW,JCW,IT1)
      AY = AY - AT1
      IF (NCR.GT.0) CALL DRSTRG(CHR,AX,AY,IC,LW,JCW,NCR)
      CALL MAPWIN(ZERO,ANX,ZERO,ANY,MNX,LNX,MNY,LNY)
C SPECIAL CASE OF VERTICAL AND HORIZONTAL LINES
      X(1) = XV
      Y(1) = ZERO
      X(2) = XV
      Y(2) = ANY
      CALL DRLINS (X,Y,2,IC,LW)
      X(1) = ZERO
      Y(1) = YV
      X(2) = ANX
      Y(2) = YV
      CALL DRLINS (X,Y,2,IC,LW)
C DRAW CONTOURS
C FIX STARTING VALUE AND INCREMENT FOR SPATIAL COORDINATES
      X(1) = ZERO
      Y(1) = ZERO
      X(2) = ONE
      Y(2) = ONE
C FIX FIRST CONTOUR LEVEL AND CONTOUR INTERVAL
      C(2) = (FMAX - FMIN)/FLOAT(NC)
      C(1) = FMIN + .5*C(2)
      CALL CONTRU (X,Y,F,C,NX,NY,NC,LINK,NXV,LW)
      IPLOT = IPLOT + 1
      IF (IPLOT.EQ.NPLOT) IPLOT = 0
      IF (IPLOT.GT.0) GO TO 70
      CALL SGRAPH
      CALL READC(IRC)
   70 RETURN
      END
*
*****************************************************
* RASTUR -- COLOR RASTER IMAGE OF FUNCTION F
*****************************************************
*
      SUBROUTINE RASTUR (F,LABEL,NX,NY,NXV,CHR,NCR,IRC)
C RASTUR DOES A COLOR RASTER IMAGE OF THE FUNCTION F, FOR VALUES OF THE
C FIRST INDEX IN THE RANGE 1 < I < N, AND OF THE SECOND INDEX 1 < J < M.
C SEVEN COLORS ARE USED.  RED, MAGENTA, YELLOW, FOREGROUND, CYAN, GREEN,
C AND BLUE ARE USED IN GOING FROM HIGHEST TO LOWEST VALUES
C ISTYLE = (0,1) = CONTOUR PLOT (FILLS AREA,PRESERVES ASPECT RATIO)
      COMMON /PLOTCM/ IPLOT, NPLOT, IDTYPE
      CHARACTER*20 LABEL
      CHARACTER*50 LBL
      DIMENSION F(NXV,NY)
      CHARACTER*(*) CHR
   91 FORMAT (A20,5H MAX=,E10.3,5H MIN=,E10.3)
      SAVE LW,ISTYLE,ZERO
      DATA LW,ISTYLE /1,1/
      DATA ZERO /0./
      IRC = 0
      IF (NPLOT.LT.1) GO TO 70
      ANX = FLOAT(NX)
      ANY = FLOAT(NY)
C GET GRAPHICS SIZE PARAMETERS
      CALL GRPARM(2,IRX,IRY,MINX,MINY,LENX,LENY,ICH,ICW,ISLB,ISTCX,ISTCY
     1,NTCX,NTCY,IGSTYL)
C FIND SCALES FOR PLOT
      FMIN = F(1,1)
      FMAX = FMIN
      DO 20 K = 1, NY
      DO 10 J = 1, NX
      IF (F(J,K).GT.FMAX) FMAX = F(J,K)
      IF (F(J,K).LT.FMIN) FMIN = F(J,K)
   10 CONTINUE
   20 CONTINUE
C FIND LOCATION FOR PLOT
      NPL1 = SQRT(FLOAT(NPLOT-1)) + 0.0001
      NPL = NPL1 + 1
      IY = IPLOT/NPL
      IX = IPLOT - IY*NPL
      MNX = (IX*IRX + MINX)/NPL
      MNY = ((NPL1 - IY)*IRY + MINY)/NPL
      LNX = LENX/NPL
      LNY = LENY/NPL
      JSLB = ISLB/NPL
      JCH = ICH/NPL
      IF (JCH.LT.1) JCH = 1
      JCW = ICW/NPL
      IF (JCW.LT.1) JCW = 1
      IF (IPLOT.GT.0) GO TO 30
      CALL INITGR(IRX,IRY,JCH,JCW)
      GO TO 40
C DRAW GRID
   30 CALL SELFMP(IRX,IRY)
   40 IC = 7
      MNX0 = MNX
      MNY0 = MNY
      IF (ISTYLE.EQ.0) GO TO 50
      MXX = MNX + LNX
      MXY = MNY + LNY
      IF (NY.GT.NX) LNX = FLOAT(LNX)*(ANX/ANY) + .5
      IF (NX.GT.NY) LNY = FLOAT(LNY)*(ANY/ANX) + .5
      MNX = MXX - LNX
      MNY = MXY - LNY
C WRITE LABELS
   50 AT1 = FLOAT(JCH + JCH/3)
      AX = FLOAT(MNX0 - JSLB)
      AY = FLOAT(MNY0) - AT1
      IT1 = 50
      WRITE (LBL,91) LABEL, FMAX, FMIN
      CALL DRSTRG(LBL,AX,AY,IC,LW,JCW,IT1)
      AY = AY - AT1
      IF (NCR.GT.0) CALL DRSTRG(CHR,AX,AY,IC,LW,JCW,NCR)
C DRAW RASTER IMAGE
      IF (FMAX.EQ.FMIN) GO TO 60
      AX = FLOAT(LNX)
      AY = FLOAT(LNY)
      CALL MAPWIN(ZERO,AX,ZERO,AY,MNX,LNX,MNY,LNY)
      CALL RASTRU(F,FMIN,FMAX,LNX,LNY,NX,NY,NXV,LW)
   60 IPLOT = IPLOT + 1
      IF (IPLOT.EQ.NPLOT) IPLOT = 0
      IF (IPLOT.GT.0) GO TO 70
      CALL SGRAPH
      CALL READC(IRC)
   70 RETURN
      END
*
*****************************************************
* TICKS  -- DRAWS BOX, TICKS, GRIDS AND LABELS AXES
*****************************************************
*
      SUBROUTINE TICKS (XMAX,XMIN,YMAX,YMIN,MINX,MINY,LENX,LENY,ICH,ICW,
     1ISLB,ISTCX,ISTCY,NTX,NTY,IC,LW,IGSTYL)
C THIS SUBROUTINE DRAWS BOX, TICKS, GRIDS AND LABELS AXES
C ICH, ICW = CHARACTER HEIGHT AND WIDTH
C NTX, NTY = NUMBER OF INTERVALS BETWEEN TICKS IN X AND Y DIRECTION
C ISTCX, ISTCY = SIZE OF TICK IN X AND Y DIRECTION
      CHARACTER*12 LBL
   91 FORMAT (E12.5)
C DRAW BOX AROUND GRAPH AND TICKS
      XMN = FLOAT(MINX)
      XMX = FLOAT(MINX + LENX)
      YMN = FLOAT(MINY)
      YMX = FLOAT(MINY + LENY)
C DRAW BOX, GRID, AND TICKS
      IF (IGSTYL.EQ.1) CALL GRIDL (XMX,XMN,YMX,YMN,ISTCX,ISTCY,NTX,NTY,I
     1C,LW,1)
      IF (IGSTYL.EQ.2) CALL GRIDG (XMX,XMN,YMX,YMN,ICH,ICW,ISTCX,ISTCY,N
     1TX,NTY,IC,LW,2)
C LABEL AXES
      AX = XMN - FLOAT(ISLB)
      AY = YMN - FLOAT(ISTCY + ICH + ICH/3)
      AT1 = XMX - FLOAT(ISLB)
      AT2 = YMX - FLOAT(ICH/2)
      IT1 = 12
      WRITE (LBL,91) XMIN
      CALL DRSTRG(LBL,XMN,AY,IC,LW,ICW,IT1)
      WRITE (LBL,91) XMAX
      CALL DRSTRG(LBL,AT1,AY,IC,LW,ICW,IT1)
      WRITE (LBL,91) YMIN
      CALL DRSTRG(LBL,AX,YMN,IC,LW,ICW,IT1)
      WRITE (LBL,91) YMAX
      CALL DRSTRG(LBL,AX,AT2,IC,LW,ICW,IT1)
      RETURN
      END
*
*****************************************************
* BOX -- DRAWS BOX
*****************************************************
*
      SUBROUTINE BOX (MINX,MINY,LENX,LENY,ICH,ICW,IC,LW,IGSTYL)
C THIS SUBROUTINE DRAWS BOX
C ICH, ICW = CHARACTER HEIGHT AND WIDTH
      XMN = FLOAT(MINX)
      XMX = FLOAT(MINX + LENX)
      YMN = FLOAT(MINY)
      YMX = FLOAT(MINY + LENY)
C DRAW BOX AROUND GRAPH
      IF (IGSTYL.EQ.1) CALL BOXL (XMX,XMN,YMX,YMN,IC,LW)
      IF (IGSTYL.EQ.2) CALL BOXG (XMX,XMN,YMX,YMN,ICH,ICW,IC,LW)
      RETURN
      END
*
*****************************************************
* CONTRU -- CONTOUR PLOT -- LOWER LEVEL ROUTINE
*****************************************************
*
      SUBROUTINE CONTRU(XA,YA,Z,C,M,N,L,LINK,MV,LWTYPE)
C SUBROUTINE WRITTEN BY ART ROSS, MODIFIED BY AEINT DE BOER FOR FORT77
C BUG FIXED BY VIKTOR DECYK.
C SEVEN COLORS ARE USED TO PLOT THE CONTOURS.  RED, MAGENTA, YELLOW,
C FOREGROUND, CYAN, GREEN, AND BLUE ARE USED IN GOING FROM HIGHEST TO
C LOWEST VALUES.
C USING RASTER GRAPHICS
      DIMENSION ICOLOR(7)
      LOGICAL*1 F1,F2,LINK
      DIMENSION C(2),LIM(2),LINK(2,M,N),XA(2),YA(2),Z(MV,N)
      EQUIVALENCE(MM1,LIM(1)),(NM1,LIM(2))
      SAVE ICOLOR
C COLORS ARE: BLUE,GREEN,CYAN,FOREGROUND,YELLOW,MAGENTA,RED
      DATA ICOLOR /1,2,3,7,6,5,4/
      MM1 = M-1
      NM1 = N-1
      AC = 7./FLOAT(L)
      DO 199 LEV=1,L
      CLEV = C(1)+(LEV-1)*C(2)
      IC = AC*(FLOAT(LEV) - .5)
      ICTYPE = ICOLOR(IC+1) + 8
C
C   MARK HORIZONTAL LINKS CROSSED BY CONTOUR.
C
      DO 10 J=1,N
      F1 = Z(1,J).GT.CLEV
      DO 10 I=1,MM1
      F2 = Z(I+1,J).GT.CLEV
      LINK(1,I,J) = F1.AND..NOT.F2.OR..NOT.F1.AND.F2
10    F1 = F2
C
C   MARK VERTICAL LINKS CROSSED BY CONTOUR.
C
      DO 20 I=1,M
      F1 = Z(I,1).GT.CLEV
      DO 20 J=1,NM1
      F2 = Z(I,J+1).GT.CLEV
      LINK(2,I,J) = F1.AND..NOT.F2.OR..NOT.F1.AND.F2
20    F1 = F2
C
C   FIRST DRAW ALL CONTOURS INTERSECTING EDGES.
C
      LX = 0
      LY = +1
      I = 1
      J = 1
      ASSIGN 101 TO IFOLLW
      DO 103 IDIR=1,4
      LIMIT = LIM(1+IABS(LX))
      LNKDIR = 1+IABS(LX)
      DO 102 K=1,LIMIT
      IF(.NOT.LINK(LNKDIR,I,J)) GO TO 101
      LINK(LNKDIR,I,J) = .FALSE.
      GO TO 1501
C
C   FOLLOWING IS INTERNAL SUBROUTINE TO FOLLOW A CONTOUR, CALLING
C THE APPROPRIATE GRAPHICS ROUTINES AS IT GOES.  IT STARTS AT THE LINK
C DESIGNATED BY I AND J, CROSSING IT IN THE DIRECTION INDICATED BY
C LX AND LY.
1501  II = I
      JJ = J
      KX = LX
      KY = LY
      ASSIGN 1502 TO IXY
      GO TO 1601
1502  CALL DRAWG(' ',XX,YY,ICTYPE,LWTYPE,0,1)
      ASSIGN 1503 TO IXY
      GO TO 1504
1503  CALL DRAWG(' ',XX,YY,ICTYPE,LWTYPE,0,2)
      LINK(1+IABS(KX),II,JJ) = .FALSE.
C
C   IL AND JL ARE INDICES OF LL CORNER OF CELL WE ARE ENTERING
C
1504  IL = II+(KX-1)/2
      JL = JJ+(KY-1)/2
      IF(IL.LT.1.OR.IL.GE.M.OR.JL.LT.1.OR.JL.GE.N)
     &  GO TO IFOLLW,(101)
      KX = -KX
      KY = -KY
      DO 1505 ICT=1,3
      ITEMP = KX
      KX = -KY
      KY = +ITEMP
      II = IL+(KX+1)/2
      JJ = JL+(KY+1)/2
      IF(LINK(1+IABS(KX),II,JJ)) GO TO 1601
1505  CONTINUE
      GO TO IFOLLW,(101)
C
C   FOLLOWING IS INTERNAL SUBROUTINE TO COMPUTE X AND Y
C COORDINATES OF THE POINT WHERE CONTOUR CROSSES THE LINK WHOSE
C INDICES ARE II AND JJ, WITH THE CROSSING IN THE DIRECTION DESIGNATED
C BY KX AND KY.  THE LINK MARK ENTRY IS NOT CLEARED.
C
1601  XX = XA(1)+(II-1)*XA(2)
      YY = YA(1)+(JJ-1)*YA(2)
      ZZ = Z(II,JJ)
      IF(KX) 1620,1610,1620
1610  X2 = XA(1)+II*XA(2)
      Z2 = Z(II+1,JJ)
      IF (Z2.NE.ZZ) XX = XX+(X2-XX)*((CLEV-ZZ)/(Z2-ZZ))
      GO TO IXY,(1502,1503)
1620  Y2 = YA(1)+JJ*YA(2)
      Z2 = Z(II,JJ+1)
      IF (Z2.NE.ZZ) YY = YY+(Y2-YY)*((CLEV-ZZ)/(Z2-ZZ))
      GO TO IXY,(1502,1503)
C
C     END OF INTERNAL SUBROUTINES
101   I = I+LY
102   J = J-LX
      ITEMP = LX
      LX = -LY
103   LY = +ITEMP
C
C   NOW DO CLOSED CONTOURS, WHICH ALL MUST EXIST ONLY IN INTERIOR OF
C ARRAY.
C
      DO 132 J=1,NM1
      DO 132 I=1,MM1
      IF(.NOT.LINK(1,I,J)) GO TO 131
      ASSIGN 131 TO IFOLLW
      LX = 0
      LY = +1
      GO TO 2501
131   IF(.NOT.LINK(2,I,J)) GO TO 132
      ASSIGN 132 TO IFOLLW
      LX = +1
      LY = 0
      GO TO 2501
C
C   FOLLOWING IS INTERNAL SUBROUTINE TO FOLLOW A CONTOUR, CALLING
C THE APPROPRIATE GRAPHICS ROUTINES AS IT GOES.  IT STARTS AT THE LINK
C DESIGNATED BY I AND J, CROSSING IT IN THE DIRECTION INDICATED BY
C LX AND LY.
2501  II = I
      JJ = J
      KX = LX
      KY = LY
      ASSIGN 2502 TO IXY
      GO TO 2601
2502  CALL DRAWG(' ',XX,YY,ICTYPE,LWTYPE,0,1)
      ASSIGN 2503 TO IXY
      GO TO 2504
2503  CALL DRAWG(' ',XX,YY,ICTYPE,LWTYPE,0,2)
      LINK(1+IABS(KX),II,JJ) = .FALSE.
C
C   IL AND JL ARE INDICES OF LL CORNER OF CELL WE ARE ENTERING
C
2504  IL = II+(KX-1)/2
      JL = JJ+(KY-1)/2
      IF(IL.LT.1.OR.IL.GE.M.OR.JL.LT.1.OR.JL.GE.N)
     &  GO TO IFOLLW,(131,132)
      KX = -KX
      KY = -KY
      DO 2505 ICT=1,3
      ITEMP = KX
      KX = -KY
      KY = +ITEMP
      II = IL+(KX+1)/2
      JJ = JL+(KY+1)/2
      IF(LINK(1+IABS(KX),II,JJ)) GO TO 2601
2505  CONTINUE
      GO TO IFOLLW,(131,132)
C
C   FOLLOWING IS INTERNAL SUBROUTINE TO COMPUTE X AND Y
C COORDINATES OF THE POINT WHERE CONTOUR CROSSES THE LINK WHOSE
C INDICES ARE II AND JJ, WITH THE CROSSING IN THE DIRECTION DESIGNATED
C BY KX AND KY.  THE LINK MARK ENTRY IS NOT CLEARED.
C
2601  XX = XA(1)+(II-1)*XA(2)
      YY = YA(1)+(JJ-1)*YA(2)
      ZZ = Z(II,JJ)
      IF(KX) 2620,2610,2620
2610  X2 = XA(1)+II*XA(2)
      Z2 = Z(II+1,JJ)
      IF (Z2.NE.ZZ) XX = XX+(X2-XX)*((CLEV-ZZ)/(Z2-ZZ))
      GO TO IXY,(2502,2503)
2620  Y2 = YA(1)+JJ*YA(2)
      Z2 = Z(II,JJ+1)
      IF (Z2.NE.ZZ) YY = YY+(Y2-YY)*((CLEV-ZZ)/(Z2-ZZ))
      GO TO IXY,(2502,2503)
C
C     END OF INTERNAL SUBROUTINES
132   CONTINUE
C
C   END OF LEVEL LOOP
C
199   CONTINUE
      RETURN
      END
*
*****************************************************
* RASTRU -- INTERNAL RASTER SUBROUTINE
*****************************************************
*
      SUBROUTINE RASTRU(F,FMIN,FMAX,LX,LY,NX,NY,NXV,LWTYPE)
C THIS SUBROUTINE CONVERTS FLOATING POINT ARRAY TO COLOR RASTER IMAGE
C SEVEN COLORS ARE USED TO CODE THE IMAGES.  RED, MAGENTA, YELLOW,
C FOREGROUND, CYAN, GREEN, AND BLUE ARE USED IN GOING FROM HIGHEST TO
C LOWEST VALUES.
C USING RASTER GRAPHICS
      DIMENSION F(NXV,NY)
      DIMENSION ICOLOR(7)
      SAVE NTC,ICOLOR
      DATA NTC /8/
C COLORS ARE: BLUE,GREEN,CYAN,FOREGROUND,YELLOW,MAGENTA,RED
      DATA ICOLOR /1,2,3,7,6,5,4/
      LX1 = LX - 1
      LY1 = LY - 1
      DXG = FLOAT(NX - 1)/FLOAT(LX1)
      DYG = FLOAT(NY - 1)/FLOAT(LY1)
      AF = FLOAT(NTC - 1)/(FMAX - FMIN)
C LOOP OVER PIXELS
      DO 40 K = 1, LY1
      Y = FLOAT(K - 1)
      YT = Y*DYG + 1.
      M = YT
      DY = YT - FLOAT(M)
      DYT = 1. - DY
      DO 20 J = 1, LX1
      X = FLOAT(J - 1)
      XT = X*DXG + 1.
      N = XT
      DX = XT - FLOAT(N)
      DXT = 1. - DX
      FC = F(N,M)*DXT*DYT + F(N+1,M)*DX*DYT + F(N,M+1)*DXT*DY + F(N+1,M+
     11)*DX*DY
      ICTYPE = 0
      IF ((FC.LT.FMIN).OR.(FC.GT.FMAX)) GO TO 10
      IC = (FC - FMIN)*AF
      IF (IC.EQ.NTC) IC = NTC - 1
      ICTYPE = ICOLOR(IC+1) + 8
   10 CALL DRAWG(' ',X,Y,ICTYPE,LWTYPE,0,3)
   20 CONTINUE
      FC = F(NX,M)*DYT + F(NX,M+1)*DY
      ICTYPE = 0
      IF ((FC.LT.FMIN).OR.(FC.GT.FMAX)) GO TO 30
      IC = (FC - FMIN)*AF
      IF (IC.EQ.NTC) IC = NTC - 1
      ICTYPE = ICOLOR(IC+1) + 8
   30 CALL DRAWG(' ',X,Y,ICTYPE,LWTYPE,0,3)
   40 CONTINUE
      FC = F(NX,NY)
      ICTYPE = 0
      IF ((FC.LT.FMIN).OR.(FC.GT.FMAX)) GO TO 50
      IC = (FC - FMIN)*AF
      IF (IC.EQ.NTC) IC = NTC - 1
      ICTYPE = ICOLOR(IC+1) + 8
   50 CALL DRAWG(' ',X,Y,ICTYPE,LWTYPE,0,3)
      RETURN
      END
*
*****************************************************
* INITGR -- INITIALIZE GRAPHICS PARAMETERS
*****************************************************
*
      SUBROUTINE INITGR(IRX,IRY,ICH,ICW)
C THIS SUBROUTINE INITIALIZES GRAPHICS PARAMETERS
C DEFAULT SCALING
      CALL SELFMP(IRX,IRY)
C DEFAULT CHARACTER SIZE
      CALL CHARSZ(ICH,ICW)
C SET CURSOR TO ZERO
      CALL DRAWG(' ',0.,0.,0,0,0,1)
C CLEAR IMAGE
      CALL DRAWG(' ',0.,0.,0,0,0,0)
      RETURN
      END
      SUBROUTINE DRLINS (X,Y,N,IC,LWTYPE)
C THIS SUBROUTINE DRAWS LINES
      DIMENSION X(N), Y(N)
      CALL DRAWG(' ',X(1),Y(1),IC,LWTYPE,0,1)
      DO 10 J = 1, N
      CALL DRAWG(' ',X(J),Y(J),IC,LWTYPE,0,2)
   10 CONTINUE
      RETURN
      END
*
*****************************************************
* DRPNTS -- DRAW POINTS
*****************************************************
*
      SUBROUTINE DRPNTS (X,Y,N,IC,LWTYPE)
C THIS SUBROUTINE DRAWS POINTS
      DIMENSION X(N), Y(N)
      DO 10 J = 1, N
      CALL DRAWG(' ',X(J),Y(J),IC,LWTYPE,0,3)
   10 CONTINUE
      RETURN
      END
*
*****************************************************
* DRSHLS -- DRAW DASHED LINES
*****************************************************
*
      SUBROUTINE DRSHLS (X,Y,N,IC,LWTYPE,L)
C THIS SUBROUTINE DRAWS DASHED LINES
      DIMENSION X(N), Y(N)
      CALL DRAWG(' ',X(1),Y(1),IC,LWTYPE,L,1)
      DO 10 J = 1, N
      CALL DRAWG(' ',X(J),Y(J),IC,LWTYPE,L,4)
   10 CONTINUE
      RETURN
      END
*
*****************************************************
* DRSTRG -- DRAW CHARACTER STRING
*****************************************************
*
      SUBROUTINE DRSTRG(CHR,AX,AY,IC,LWTYPE,ICW,NCR)
C THIS SUBROUTINE DRAWS CHARACTER STRING OF LENGTH NCR
      CHARACTER*(*) CHR
      DX = FLOAT(ICW)
      IF (NCR.EQ.0) GO TO 20
      DO 10 I = 1, NCR
      AT1 = AX + DX*FLOAT(I - 1)
      CALL DRAWG(CHR(I:I),AT1,AY,IC,LWTYPE,0,5)
   10 CONTINUE
   20 RETURN
      END
*
*****************************************************
* GRIDL -- DRAW TICKS AND/OR GRIDS
*****************************************************
*
      SUBROUTINE GRIDL (XMX,XMN,YMX,YMN,ISTCX,ISTCY,NTX,NTY,IC,LWTYPE,IS
     1TYLE)
C THIS SUBROUTINE DRAWS TICKS AND/OR GRIDS ON GRAPH WITH LINES
C ISTYLE = (0,1,2) = DRAW (BOX,TICKS,TICKS AND GRID)
C DRAW TICKS IN Y DIRECTION
      IT1 = NTY + 1
      STX = FLOAT(ISTCX)
      AT1 = XMN - STX
      AT2 = XMX + STX
      DYT = (YMX - YMN)/FLOAT(NTY)
      DO 20 J = 1, IT1
      AY = DYT*FLOAT(J - 1) + YMN
      CALL DRAWG(' ',AT1,AY,IC,LWTYPE,0,1)
      IF ((J.EQ.1).OR.(J.EQ.IT1)) GO TO 10
      IF (ISTYLE.EQ.0) GO TO 20
      CALL DRAWG(' ',XMN,AY,IC,LWTYPE,0,4)
      IF (ISTYLE.EQ.1) CALL DRAWG(' ',XMX,AY,IC,LWTYPE,0,1)
      IF (ISTYLE.EQ.2) CALL DRAWG(' ',XMX,AY,IC,LWTYPE,1,4)
   10 CALL DRAWG(' ',AT2,AY,IC,LWTYPE,0,4)
   20 CONTINUE
C DRAW TICKS IN X DIRECTION
      IT1 = NTX + 1
      STY = FLOAT(ISTCY)
      AT1 = YMN - STY
      AT2 = YMX + STY
      DXT = (XMX - XMN)/FLOAT(NTX)
      DO 40 J = 1, IT1
      AX = DXT*FLOAT(J - 1) + XMN
      CALL DRAWG(' ',AX,AT1,IC,LWTYPE,0,1)
      IF ((J.EQ.1).OR.(J.EQ.IT1)) GO TO 30
      IF (ISTYLE.EQ.0) GO TO 40
      CALL DRAWG(' ',AX,YMN,IC,LWTYPE,0,4)
      IF (ISTYLE.EQ.1) CALL DRAWG(' ',AX,YMX,IC,LWTYPE,0,1)
      IF (ISTYLE.EQ.2) CALL DRAWG(' ',AX,YMX,IC,LWTYPE,1,4)
   30 CALL DRAWG(' ',AX,AT2,IC,LWTYPE,0,4)
   40 CONTINUE
      RETURN
      END
*
*****************************************************
* GRIDG  -- DRAW TICKS AND/OR GRIDS
*****************************************************
*
      SUBROUTINE GRIDG (XMX,XMN,YMX,YMN,ICH,ICW,ISTCX,ISTCY,NTX,NTY,IC,L
     1WTYPE,ISTYLE)
C THIS SUBROUTINE DRAWS TICKS AND/OR GRIDS ON GRAPH WITH CHARACTERS
C ISTYLE = (0,1,2) = DRAW (BOX,TICKS,TICKS AND GRID)
      DX = FLOAT(ICW)
      DY = FLOAT(ICH)
      DXH = FLOAT(ICW/2)
      DYH = FLOAT(ICH/2)
      DXT = (XMX - XMN)/FLOAT(NTX)
      DYT = (YMX - YMN)/FLOAT(NTY)
      NCRX = (XMX - XMN)/DX + 1.5
      NCRY = (YMX - YMN)/DY + 1.5
      DX = (XMX - XMN)/FLOAT(NCRX - 1)
      DY = (YMX - YMN)/FLOAT(NCRY - 1)
C DRAW TICKS IN Y DIRECTION
      IT1 = NTY + 1
      STX = FLOAT(ISTCX)
      AT1 = XMN - STX
      AT2 = XMX + STX
      NCRT = STX/DX + .5
      DO 70 J = 1, IT1
      IN = 1
      IF ((J.EQ.1).OR.(J.EQ.IT1)) IN = 0
      IY = DYT*FLOAT(J - 1)/DY + .5
      AY = DY*FLOAT(IY) + YMN - DYH
C EXTERIOR LEFT TICKS
      IF ((IN.EQ.1).AND.(ISTYLE.EQ.0)) GO TO 20
      DO 10 I = 1, NCRT
      AT3 = AT1 + DX*FLOAT(I - 1) - DXH
      CALL DRAWG('-',AT3,AY,IC,LWTYPE,0,5)
   10 CONTINUE
C GRIDS AND INTERIOR TICKS
   20 AX = DXT
      IX = AX/DX + .5
      DO 50 I = 1, NCRX
      I1 = I - 1
      AT3 = XMN + DX*FLOAT(I1) - DXH
      IF ((I.EQ.1).OR.(I.EQ.NCRX)) GO TO 40
      IF ((IN.EQ.1).AND.(ISTYLE.LE.1)) GO TO 50
      IF ((ISTYLE.GT.0).AND.(I1.EQ.IX)) GO TO 30
      CALL DRAWG('-',AT3,AY,IC,LWTYPE,0,5)
      GO TO 50
   30 AX = AX + DXT
      IX = AX/DX + .5
   40 CALL DRAWG('+',AT3,AY,IC,LWTYPE,0,5)
   50 CONTINUE
C EXTERIOR RIGHT TICKS
      IF ((IN.EQ.1).AND.(ISTYLE.EQ.0)) GO TO 70
      DO 60 I = 1, NCRT
      AT3 = XMX + DX*FLOAT(I) - DXH
      CALL DRAWG('-',AT3,AY,IC,LWTYPE,0,5)
   60 CONTINUE
   70 CONTINUE
C DRAW TICKS IN X DIRECTION
      IT1 = NTX + 1
      STY = FLOAT(ISTCY)
      AT1 = YMN - STY
      AT2 = YMX + STY
      NCRT = STY/DY + .5
      DO 130 J = 1, IT1
      IN = 1
      IF ((J.EQ.1).OR.(J.EQ.IT1)) IN = 0
      IX = DXT*FLOAT(J - 1)/DX + .5
      AX = DX*FLOAT(IX) + XMN - DXH
C EXTERIOR BOTTOM TICKS
      IF ((IN.EQ.1).AND.(ISTYLE.EQ.0)) GO TO 90
      DO 80 I = 1, NCRT
      AT3 = AT1 + DY*FLOAT(I - 1) - DYH
      CALL DRAWG('|',AX,AT3,IC,LWTYPE,0,5)
   80 CONTINUE
C GRIDS AND INTERIOR TICKS
   90 AY = DYT
      IY = AY/DY + .5
      DO 110 I = 1, NCRY
      I1 = I - 1
      AT3 = YMN + DY*FLOAT(I1) - DYH
      IF ((I.EQ.1).OR.(I.EQ.NCRY)) GO TO 110
      IF ((IN.EQ.1).AND.(ISTYLE.LE.1)) GO TO 110
      IF ((ISTYLE.GT.0).AND.(I1.EQ.IY)) GO TO 100
      CALL DRAWG('|',AX,AT3,IC,LWTYPE,0,5)
      GO TO 110
  100 AY = AY + DYT
      IY = AY/DY + .5
  110 CONTINUE
C EXTERIOR TOP TICKS
      IF ((IN.EQ.1).AND.(ISTYLE.EQ.0)) GO TO 130
      DO 120 I = 1, NCRT
      AT3 = YMX + DY*FLOAT(I) - DYH
      CALL DRAWG('|',AX,AT3,IC,LWTYPE,0,5)
  120 CONTINUE
  130 CONTINUE
      RETURN
      END
*
*****************************************************
* BOXL   -- DRAW BOX AROUND GRAPH WITH LINES
*****************************************************
*
      SUBROUTINE BOXL (XMX,XMN,YMX,YMN,IC,LWTYPE)
C THIS SUBROUTINE DRAWS BOX AROUND GRAPH WITH LINES
      CALL DRAWG(' ',XMN,YMN,IC,LWTYPE,0,1)
      CALL DRAWG(' ',XMX,YMN,IC,LWTYPE,0,2)
      CALL DRAWG(' ',XMX,YMX,IC,LWTYPE,0,2)
      CALL DRAWG(' ',XMN,YMX,IC,LWTYPE,0,2)
      CALL DRAWG(' ',XMN,YMN,IC,LWTYPE,0,2)
      RETURN
      END
*
*****************************************************
* BOXG   -- DRAW BOX AROUND GRAPH WITH CHARACTERS
*****************************************************
*
      SUBROUTINE BOXG (XMX,XMN,YMX,YMN,ICH,ICW,IC,LWTYPE)
C THIS SUBROUTINE DRAWS BOX AROUND GRAPH WITH CHARACTERS
      DX = FLOAT(ICW)
      DY = FLOAT(ICH)
      DXH = FLOAT(ICW/2)
      DYH = FLOAT(ICH/2)
      NCRX = (XMX - XMN)/DX + 1.5
      NCRY = (YMX - YMN)/DY + 1.5
      DX = (XMX - XMN)/FLOAT(NCRX - 1)
      DY = (YMX - YMN)/FLOAT(NCRY - 1)
C DRAWS LINES IN X DIRECTION
      AT2 = YMN - DYH
      AT3 = YMX - DYH
      DO 20 I = 1, NCRX
      AT1 = XMN + DX*FLOAT(I - 1) - DXH
      IF ((I.EQ.1).OR.(I.EQ.NCRX)) GO TO 10
      CALL DRAWG('-',AT1,AT2,IC,LWTYPE,0,5)
      CALL DRAWG('-',AT1,AT3,IC,LWTYPE,0,5)
      GO TO 20
   10 CALL DRAWG('+',AT1,AT2,IC,LWTYPE,0,5)
      CALL DRAWG('+',AT1,AT3,IC,LWTYPE,0,5)
   20 CONTINUE
C DRAWS LINES IN Y DIRECTION
      AT2 = XMN - DXH
      AT3 = XMX - DXH
      DO 30 I = 1, NCRY
      IF ((I.EQ.1).OR.(I.EQ.NCRY)) GO TO 30
      AT1 = YMN + DY*FLOAT(I - 1) - DYH
      CALL DRAWG('|',AT2,AT1,IC,LWTYPE,0,5)
      CALL DRAWG('|',AT3,AT1,IC,LWTYPE,0,5)
   30 CONTINUE
      RETURN
      END
*
*****************************************************
* SGRAPH -- WRITES OUT PLOT TO DEVICE
*****************************************************
*
      SUBROUTINE SGRAPH
C THIS SUBROUTINE WRITES OUT PLOT TO DEVICE
      CALL DRAWG(' ',0.,0.,0,0,0,6)
      RETURN
      END
*
*****************************************************
* READC
*****************************************************
*
      SUBROUTINE READC(IRC)
      CHARACTER*1 C
      CHARACTER*8 CX
      CHARACTER*37 CHR1
      CHARACTER*13 CHR2
   91 FORMAT (1X,A37)
   92 FORMAT (1X,A13)
      DATA CHR1 /' Q=QUIT, S=SAVE, M=MODIFY, R=REVERSE '/
      DATA CHR2 /' CR=CONTINUE '/
      IRC = 0
   10 CALL CINPUT(CX)
      C = CX(1:1)
      IF (C.NE.'?') GO TO 20
      CALL CLEAR
      WRITE (6,91) CHR1
      WRITE (6,92) CHR2
      GO TO 10
   20 IF ((C.EQ.'Q').OR.(C.EQ.'q')) IRC = 1
      IF ((C.EQ.'S').OR.(C.EQ.'s')) IRC = 2
      IF ((C.EQ.'M').OR.(C.EQ.'m')) IRC = 3
      IF ((C.EQ.'R').OR.(C.EQ.'r')) IRC = 4
      RETURN
      END
*
*****************************************************
* WPARAM --
*****************************************************
*
      SUBROUTINE WPARAM (RUNID,INDX,MOVION,NPX,NPXB,T0,TEND,DT,AX,VTX,ED
     1GE,RMASS,RTEMP,VDX,VTDX,IBCS,BVL,BVR,NUSTRT,NTR,NTI,NPP,NTW,NTT,NP
     2ROBT,NDIST,NTP,MODES,NDP,NTV,NMV,NXB,NPRS,NPRO,NDV,NTD,MODED,NDD,A
     3NLE,ANSE,AMODEN,FREQN,ANLI,ANSI,AMODEX,FREQ,TRMP,TOFF,QME,QMI,QMB,
     4QTEST,VTEST,X0,CI,IRC)
      CHARACTER*128 CHR
      CHARACTER*8 RUNID
  901 FORMAT (8H RUNID= ,A8)
  902 FORMAT (6H INDX=,I3,5H NPX=,I8,5H VTX=,F10.7,4H CI=,F8.4,5H QME=,F
     18.4)
  903 FORMAT (6H IBCS=,I2,5H BVL=,F14.7,5H BVR=,F14.7)
  904 FORMAT (4H T0=,F8.1,6H TEND=,F8.1,4H DT=,F8.5)
  905 FORMAT (8H MOVION=,I2,7H RMASS=,F14.7,7H RTEMP=,F14.7,5H QMI=,F8.4
     1)
  906 FORMAT (6H NPXB=,I8,5H VDX=,F10.7,6H VTDX=,F10.7,5H QMB=,F8.4)
  907 FORMAT (5H NPP=,I6,5H NTI=,I4,4H AX=,F8.5,6H EDGE=,F8.5)
  908 FORMAT (8H NUSTRT=,I2,5H NTR=,I6)
  909 FORMAT (5H NTW=,I6)
  910 FORMAT (5H NTT=,I6,8H NPROBT=,I8,7H NDIST=,I5)
  911 FORMAT (5H NTP=,I6,7H MODES=,I5,5H NDP=,I6)
  912 FORMAT (5H NTV=,I6,5H NMV=,I5,5H NXB=,I4,6H NPRS=,I3,6H NPRO=,I3,5
     1H NDV=,I6)
  913 FORMAT (5H NTD=,I6,7H MODED=,I5,5H NDD=,I6)
  914 FORMAT (6H ANLE=,F14.7,6H ANSE=,F14.7)
  915 FORMAT (8H AMODEN=,F8.3,7H FREQN=,F14.7)
  916 FORMAT (6H ANLI=,F14.7,6H ANSI=,F14.7)
  917 FORMAT (8H AMODEX=,F8.3,6H FREQ=,F14.7,6H TRMP=,F8.1,6H TOFF=,F8.1
     1)
  918 FORMAT (7H QTEST=,F11.5,7H VTEST=,F5.2,4H X0=,F8.1)
      SAVE LW,IC
      DATA LW,IC /1,7/
      IRC = 0
      CALL GRPARM(3,IRX,IRY,MINX,MINY,LENX,LENY,ICH,ICW,ISLB,ISTCX,ISTCY
     1,NTCX,NTCY,IGSTYL)
      CALL INITGR(IRX,IRY,ICH,ICW)
      AT1 = FLOAT(ICH + ICH/3)
      AX = FLOAT(MINX)
      WRITE (CHR,901) RUNID
      AY = FLOAT(MINY + LENY) - AT1
      CALL DRSTRG(CHR,AX,AY,IC,LW,ICW,16)
      WRITE (CHR,902) INDX, NPX, VTX, CI, QME
      AY = AY - AT1
      CALL DRSTRG(CHR,AX,AY,IC,LW,ICW,62)
      WRITE (CHR,903) IBCS, BVL, BVR
      AY = AY - AT1
      CALL DRSTRG(CHR,AX,AY,IC,LW,ICW,46)
      WRITE (CHR,904) T0, TEND, DT
      AY = AY - AT1
      CALL DRSTRG(CHR,AX,AY,IC,LW,ICW,38)
      WRITE (CHR,905) MOVION, RMASS, RTEMP, QMI
      AY = AY - AT1
      CALL DRSTRG(CHR,AX,AY,IC,LW,ICW,65)
      WRITE (CHR,906) NPXB, VDX, VTDX, QMB
      AY = AY - AT1
      CALL DRSTRG(CHR,AX,AY,IC,LW,ICW,58)
      WRITE (CHR,907) NPP, NTI, AX, EDGE
      AY = AY - AT1
      CALL DRSTRG(CHR,AX,AY,IC,LW,ICW,46)
      WRITE (CHR,908) NUSTRT, NTR
      AY = AY - AT1
      CALL DRSTRG(CHR,AX,AY,IC,LW,ICW,21)
      WRITE (CHR,909) NTW
      AY = AY - AT1
      CALL DRSTRG(CHR,AX,AY,IC,LW,ICW,11)
      WRITE (CHR,910) NTT, NPROBT, NDIST
      AY = AY - AT1
      CALL DRSTRG(CHR,AX,AY,IC,LW,ICW,39)
      WRITE (CHR,911) NTP, MODES, NDP
      AY = AY - AT1
      CALL DRSTRG(CHR,AX,AY,IC,LW,ICW,34)
      WRITE (CHR,912) NTV, NMV, NXB, NPRS, NPRO, NDV
      AY = AY - AT1
      CALL DRSTRG(CHR,AX,AY,IC,LW,ICW,60)
      WRITE (CHR,913) NTD, MODED, NDD
      AY = AY - AT1
      CALL DRSTRG(CHR,AX,AY,IC,LW,ICW,34)
      WRITE (CHR,914) ANLE, ANSE
      AY = AY - AT1
      CALL DRSTRG(CHR,AX,AY,IC,LW,ICW,40)
      WRITE (CHR,915) AMODEN, FREQN
      AY = AY - AT1
      CALL DRSTRG(CHR,AX,AY,IC,LW,ICW,37)
      WRITE (CHR,916) ANLI, ANSI
      AY = AY - AT1
      CALL DRSTRG(CHR,AX,AY,IC,LW,ICW,40)
      WRITE (CHR,917) AMODEX, FREQ, TRMP, TOFF
      AY = AY - AT1
      CALL DRSTRG(CHR,AX,AY,IC,LW,ICW,64)
      WRITE (CHR,918) QTEST, VTEST, X0
      AY = AY - AT1
      CALL DRSTRG(CHR,AX,AY,IC,LW,ICW,42)
      CALL SGRAPH
      CALL READC(IRC)
      RETURN
      END
*
*****************************************************
* WPCORR
*****************************************************
*
      SUBROUTINE WPCORR (RUNID,INDX,NTP,MODES,IBCS,T0,TEND,DT,CENG,LTS,I
     1TS,NTS,KMIN,KMAX,NTD,NTC,WMIN,WMAX,DW,IRC)
      CHARACTER*128 CHR
      CHARACTER*8 RUNID
  900 FORMAT (39H SPECTRUM ANALYSIS FOR 1D PERIODIC DATA)
  901 FORMAT (8H RUNID= ,A8,6H INDX=,I3,5H NTP=,I6,7H MODES=,I5,6H IBCS=
     1,I2)
  902 FORMAT (4H T0=,F8.1,6H TEND=,F8.1,4H DT=,F8.5,6H CENG=,E14.7)
  903 FORMAT (5H LTS=,I6,5H ITS=,I6, 5H NTS=,I6)
  904 FORMAT (6H KMIN=,I6,6H KMAX=,I6)
  905 FORMAT (5H NTD=,I6,5H NTC=,I6)
  906 FORMAT (6H WMIN=,F8.4,6H WMAX=,F8.4,4H DW=,F8.4)
      SAVE LW,IC
      DATA LW,IC /1,7/
      IRC = 0
      CALL GRPARM(3,IRX,IRY,MINX,MINY,LENX,LENY,ICH,ICW,ISLB,ISTCX,ISTCY
     1,NTCX,NTCY,IGSTYL)
      CALL INITGR(IRX,IRY,ICH,ICW)
      AT1 = FLOAT(ICH + ICH/3)
      AX = FLOAT(MINX)
      WRITE (CHR,900)
      AY = FLOAT(MINY + LENY) - 2.*AT1
      CALL DRSTRG(CHR,AX,AY,IC,LW,ICW,39)
      WRITE (CHR,901) RUNID, INDX, NTP, MODES, IBCS
      AY = AY - AT1
      CALL DRSTRG(CHR,AX,AY,IC,LW,ICW,56)
      WRITE (CHR,902) T0, TEND, DT, CENG
      AY = AY - AT1
      CALL DRSTRG(CHR,AX,AY,IC,LW,ICW,58)
      WRITE (CHR,903) LTS, ITS, NTS
      AY = AY - AT1
      CALL DRSTRG(CHR,AX,AY,IC,LW,ICW,33)
      WRITE (CHR,904) KMIN, KMAX
      AY = AY - AT1
      CALL DRSTRG(CHR,AX,AY,IC,LW,ICW,24)
      WRITE (CHR,905) NTD, NTC
      AY = AY - AT1
      CALL DRSTRG(CHR,AX,AY,IC,LW,ICW,22)
      WRITE (CHR,906) WMIN, WMAX, DW
      AY = AY - AT1
      CALL DRSTRG(CHR,AX,AY,IC,LW,ICW,40)
      CALL SGRAPH
      CALL READC(IRC)
      RETURN
      END
*
*****************************************************
* WPCLDS
*****************************************************
*
      SUBROUTINE WPCLDS (RUNID,INDX,NTP,MODES,IBCS,T0,TEND,DT,CENG,VTEST
     1,QTEST,NP,LTS,ITS,NTS,MTS,NXD,NXS,LAB,IRC)
      CHARACTER*128 CHR
      CHARACTER*8 RUNID
  900 FORMAT (38H DISPLAY SUBTRACTED DATA FOR 1D CLOUDS)
  901 FORMAT (8H RUNID= ,A8,6H INDX=,I3,5H NTP=,I6,7H MODES=,I5,6H IBCS=
     1,I2)
  902 FORMAT (4H T0=,F8.1,6H TEND=,F8.1,4H DT=,F8.5,6H CENG=,E14.7)
  903 FORMAT (9H VTEST = ,F5.2,9H QTEST = ,F11.5,6H NP = ,I8)
  904 FORMAT (5H LTS=,I6,5H ITS=,I6,5H NTS=,I6,5H MTS=,I6)
  905 FORMAT (5H LAB=,I6,5H NXD=,I6,5H NXS=,I6)
      SAVE LW,IC
      DATA LW,IC /1,7/
      IRC = 0
      CALL GRPARM(3,IRX,IRY,MINX,MINY,LENX,LENY,ICH,ICW,ISLB,ISTCX,ISTCY
     1,NTCX,NTCY,IGSTYL)
      CALL INITGR(IRX,IRY,ICH,ICW)
      AT1 = FLOAT(ICH + ICH/3)
      AX = FLOAT(MINX)
      WRITE (CHR,900)
      AY = FLOAT(MINY + LENY) - 2.*AT1
      CALL DRSTRG(CHR,AX,AY,IC,LW,ICW,38)
      WRITE (CHR,901) RUNID, INDX, NTP, MODES, IBCS
      AY = AY - AT1
      CALL DRSTRG(CHR,AX,AY,IC,LW,ICW,56)
      WRITE (CHR,902) T0, TEND, DT, CENG
      AY = AY - AT1
      CALL DRSTRG(CHR,AX,AY,IC,LW,ICW,58)
      WRITE (CHR,903) VTEST, QTEST, NP
      AY = AY - AT1
      CALL DRSTRG(CHR,AX,AY,IC,LW,ICW,48)
      WRITE (CHR,904) LTS, ITS, NTS, MTS
      AY = AY - AT1
      CALL DRSTRG(CHR,AX,AY,IC,LW,ICW,44)
      WRITE (CHR,905) LAB, NXD, NXS
      AY = AY - AT1
      CALL DRSTRG(CHR,AX,AY,IC,LW,ICW,33)
      CALL SGRAPH
      CALL READC(IRC)
      RETURN
      END
*
*****************************************************
* WPPRFL
*****************************************************
*
      SUBROUTINE WPPRFL (RUNID,INDX,MOVION,NTV,NPRS,NPRO,T0,TEND,DT,RMAS
     1S,LTS,ITS,NTS,MTS,IDS,NDS,ION,IFL,IRC)
      CHARACTER*128 CHR
      CHARACTER*8 RUNID
  900 FORMAT (37H DISPLAY SPATIAL PROFILES FOR 1D DATA)
  901 FORMAT (8H RUNID= ,A8,6H INDX=,I3,8H MOVION=,I2,5H NTV=,I6,6H NPRS
     1=,I3,6H NPRO=,I3)
  902 FORMAT (4H T0=,F8.1,6H TEND=,F8.1,4H DT=,F8.5,7H RMASS=,F14.7)
  903 FORMAT (5H LTS=,I6,5H ITS=,I6,5H NTS=,I6,5H MTS=,I6)
  904 FORMAT (5H IDS=,I6,5H NDS=,I6,5H ION=,I6,5H IFL=,I6)
      SAVE LW,IC
      DATA LW,IC /1,7/
      IRC = 0
      CALL GRPARM(3,IRX,IRY,MINX,MINY,LENX,LENY,ICH,ICW,ISLB,ISTCX,ISTCY
     1,NTCX,NTCY,IGSTYL)
      CALL INITGR(IRX,IRY,ICH,ICW)
      AT1 = FLOAT(ICH + ICH/3)
      AX = FLOAT(MINX)
      WRITE (CHR,900)
      AY = FLOAT(MINY + LENY) - 2.*AT1
      CALL DRSTRG(CHR,AX,AY,IC,LW,ICW,37)
      WRITE (CHR,901) RUNID, INDX, MOVION, NTV, NPRS, NPRO
      AY = AY - AT1
      CALL DRSTRG(CHR,AX,AY,IC,LW,ICW,64)
      WRITE (CHR,902) T0, TEND, DT, RMASS
      AY = AY - AT1
      CALL DRSTRG(CHR,AX,AY,IC,LW,ICW,59)
      WRITE (CHR,903) LTS, ITS, NTS, MTS
      AY = AY - AT1
      CALL DRSTRG(CHR,AX,AY,IC,LW,ICW,44)
      WRITE (CHR,904) IDS, NDS, ION, IFL
      AY = AY - AT1
      CALL DRSTRG(CHR,AX,AY,IC,LW,ICW,44)
      CALL SGRAPH
      CALL READC(IRC)
      RETURN
      END
*
*****************************************************
* GCLOSE
*****************************************************
*
      SUBROUTINE GCLOSE
C THIS SUBROUTINE CLOSES GRAPHICS LIBRARY
      COMMON /PLOTCM/ IPLOT, NPLOT, IDTYPE
      IF (IPLOT.EQ.0) GO TO 10
      CALL SGRAPH
      CALL READC(IRC)
   10 CALL QUITG
      RETURN
      END
*
*****************************************************
* VERSTC -- TEKTRONIX TO RASTER CONVERSION
*****************************************************
*
C TEKTRONIX EMULATOR LIBRARY FOR RASTER PLOTS
      SUBROUTINE VERSTC (IA,IC,IRB)
C THIS PROGRAM SENDS TEKTRONIX 4012 INFORMATION TO RASTER FILE
C WRITTEN FOR THE IBM 3090VF - VIKTOR K. DECYK, UCLA
      COMMON /DEVICE/ ID, ICFLG
      CHARACTER*1 LBL
      DIMENSION IA(4)
      DIMENSION LT(8), ICOLOR(8), LW(3), IATE(128)
      SAVE IG,NF,ND,IRD,IX,IY,IF,LS,LT,ICOLOR,LW,IATE
      SAVE XS,YS,ICX,ICY,IHY,IBY,LTYPE,ICTYPE,LWTYPE
      SAVE LX,LY,IYH,IYL,IXH,IXL
   93 FORMAT (1H1,I6,17H FRAME(S) PLOTTED)
      DATA IG,NF,ND,IP,IRD,IZ,IF,LS /0,0,0,0,0,0,0,0/
C LINE STYLE TABLE
      DATA LT /0,1,2,3,4,0,0,0/
C COLORS TABLE: FOREGROUND,BLUE,RED,YELLOW,CYAN,MAGENTA,GREEN,FOREGROUND
      DATA ICOLOR /7,1,4,6,3,5,2,7/
C LINE WIDTH TABLE
      DATA LW /1,2,1/
C EBCDIC CODE FOR ASCII 124 IS NON-STANDARD
      DATA IATE /0,1,2,3,55,45,46,47,22,5,37,11,12,13,14,15,16,17,18,19,
     160,61,50,38,24,25,63,39,34,29,53,31,64,90,127,123,91,108,80,125,77
     2,93,92,78,107,96,75,97,240,241,242,243,244,245,246,247,248,249,122
     3,94,76,126,110,111,124,193,194,195,196,197,198,199,200,201,209,210
     4,211,212,213,214,215,216,217,226,227,228,229,230,231,232,233,173,2
     524,189,95,109,121,129,130,131,132,133,134,135,136,137,145,146,147,
     6148,149,150,151,152,153,162,163,164,165,166,167,168,169,192,79,208
     7,161,7/
      IRB = 0
      IF (IC.LT.1) GO TO 300
      K = 1
      IF (IRD.EQ.0) GO TO 290
      GO TO (10,70,220,230,250), IRD
C READ CHARACTER
   10 IRD = 1
      IF (K.GT.IC) GO TO 310
      LC = IA(K)
      K = K + 1
   20 IF (LC.LT.32) GO TO 30
      IF (IG.GT.0) GO TO 200
C ALPHA MODE
C     LBL = CHAR(LC)
      LBL = CHAR(IATE(LC+1))
C PLOT CHARACTER
      AX = FLOAT(IX)
      AY = FLOAT(IY)
      CALL DRAWG(LBL,AX,AY,ICTYPE,LWTYPE,0,5)
      GO TO 150
C CONTROL CHARACTERS
   30 IF (LC.EQ.29) GO TO 40
      IF (LC.EQ.13) GO TO 50
      IF (LC.EQ.31) GO TO 60
      IF (LC.EQ.27) GO TO 70
      IF (LC.EQ.28) GO TO 140
      IF (LC.EQ.9) GO TO 150
      IF (LC.EQ.10) GO TO 160
      IF (LC.EQ.8) GO TO 170
      IF (LC.EQ.11) GO TO 180
      IF (LC.EQ.7) GO TO 190
C UNKNOWN CONTROL CHARACTER
      GO TO 10
C SET GRAPH MODE (DARK VECTOR)
   40 IG = 2
      GO TO 10
C CARRIAGE RETURN
   50 IX = IZ
C SET ALPHA MODE
   60 IG = 0
      GO TO 10
C ESCAPE SEQUENCE
   70 IRD = 2
      IF (K.GT.IC) GO TO 310
      LC = IA(K)
      K = K + 1
      IF (LC.EQ.12) GO TO 80
      IF (LC.EQ.23) GO TO 90
      IF (LC.EQ.56) GO TO 100
      IF (LC.EQ.57) GO TO 110
      IF (LC.EQ.58) GO TO 120
      IF (LC.EQ.59) GO TO 130
      IF ((LC.GE.96).AND.(LC.LE.119)) GO TO 135
C UNSUPPORTED ESCAPE SEQUENCE
      GO TO 10
C NEW FRAME
   80 NF = NF + 1
      IF (NF.LT.ND) GO TO 89
      IF (NF.EQ.ND) GO TO 85
      IF (ID.EQ.1) CALL LABELF(NF)
      CALL DRAWG(' ',0.,0.,0,0,0,6)
      ND = ND + 1
      IF (ID.GT.1) GO TO 85
      CALL READN(IRC,NRC)
      IF (IRC.NE.1) GO TO 83
      CALL QUITG
      WRITE (6,93) NF
      STOP 1
   83 IF (IRC.EQ.2) ND = NRC - 1
      IF (NF.LT.ND) GO TO 89
      IF (NF.EQ.ND) GO TO 85
      NF = 0
      IRD = 1
      IRB = 1
      IF (NF.LT.ND) GO TO 84
      CALL CHARSZ(ICY,ICX)
      CALL DRAWG(' ',0.,0.,0,0,0,0)
      LTYPE = LT(1)
      ICTYPE = ICOLOR(1)
      LWTYPE = LW(1)
   84 IG = 0
      LS = 0
      IF = 0
      IX = IZ
      IY = IHY
      GO TO 310
   85 CALL DRAWG(' ',0.,0.,0,0,0,0)
      LTYPE = LT(1)
      ICTYPE = ICOLOR(1)
      LWTYPE = LW(1)
   89 IG = 0
      LS = 0
      IF = 0
      IX = IZ
      IY = IHY
      GO TO 10
C MAKE HARDCOPY
   90 IF (IP.GT.0) GO TO 10
      GO TO 10
C LARGE CHARACTERS
  100 ICX = 14.*XS + .5
      ICY = 22.*YS + .5
      CALL CHARSZ(ICY,ICX)
      GO TO 10
C MEDIUM-LARGE CHARACTERS
  110 ICX = 13.*XS + .5
      ICY = 21.*YS + .5
      CALL CHARSZ(ICY,ICX)
      GO TO 10
C MEDIUM-SMALL CHARACTERS
  120 ICX = 9.*XS + .5
      ICY = 13.*YS + .5
      CALL CHARSZ(ICY,ICX)
      GO TO 10
C SMALL CHARACTERS
  130 ICX = 8.*XS + .5
      ICY = 12.*YS + .5
      CALL CHARSZ(ICY,ICX)
      GO TO 10
C SET LINE STYLE AND FOCUS
  135 LS = LC - 96
      IF = LS/8
      LS = LS - IF*8
      LTYPE = LT(1)
      IF (ICFLG.NE.1) LTYPE = LT(LS+1)
      ICTYPE = ICOLOR(1)
      IF (ICFLG.NE.0) ICTYPE = ICOLOR(LS+1)
      LWTYPE = LW(IF+1)
      GO TO 10
C POINT PLOTTING MODE
  140 IG = 3
      GO TO 10
C TAB
  150 IX = IX + ICX
      IF (IX.LT.LX) GO TO 10
      IX = IX - LX
C LINE FEED
  160 IY = IY - ICY
      IF (IY.LT.IZ) IY = IHY
      GO TO 10
C BACKSPACE
  170 IX = IX - ICX
      IF (IX.GE.IZ) GO TO 10
      IX = IX + LX
C VERTICAL TAB
  180 IY = IY + ICY
      IF (IY.GE.LY) IY = IBY
      GO TO 10
C BELL (SET VECTOR TO DRAW)
  190 IF (IG.EQ.2) IG = 1
      GO TO 10
C GRAPH MODE
C DECODE ADDRESS
  200 IF (LC.LT.64) GO TO 240
      IF (LC.LT.96) GO TO 260
  210 IYL = LC
  220 IRD = 3
      IF (K.GT.IC) GO TO 310
      LC = IA(K)
      K = K + 1
      IF (LC.GE.96) GO TO 210
      IF (LC.GE.64) GO TO 260
      IF (LC.LT.32) GO TO 220
      IXH = LC
  230 IRD = 4
      IF (K.GT.IC) GO TO 310
      LC = IA(K)
      K = K + 1
      IF (LC.GE.32) GO TO 260
      GO TO 230
  240 IYH = LC
  250 IRD = 5
      IF (K.GT.IC) GO TO 310
      LC = IA(K)
      K = K + 1
      IF (LC.GE.96) GO TO 210
      IF (LC.LT.32) GO TO 250
  260 IXL = LC
C CALCULATE ADDRESS
      I = (IXH - 34)*32 + IXL
      J = (IYH - 35)*32 + IYL
      I = XS*FLOAT(I) + .5
      J = YS*FLOAT(J) + .5
C PERFORM CLIPPING
      IF (I.LT.IZ) I = IZ
      IF (I.GE.LX) I = LX - 1
      IF (J.LT.IZ) J = IZ
      IF (J.GE.LY) J = LY - 1
C CONVERT TO FLOATING POINT REPRESENTATION
      AX = FLOAT(I)
      AY = FLOAT(J)
      IF (IG.NE.1) GO TO 270
C DRAW VECTOR
      IF (ICFLG.EQ.1) CALL DRAWG(' ',AX,AY,ICTYPE,LWTYPE,0,2)
      IF (ICFLG.NE.1) CALL DRAWG(' ',AX,AY,ICTYPE,LWTYPE,LTYPE,4)
      IX = I
      IY = J
      GO TO 10
C PERFORM MOVE AND SET VECTOR TO DRAW
  270 IF (IG.EQ.3) GO TO 280
      CALL DRAWG(' ',AX,AY,ICTYPE,LWTYPE,0,1)
      IX = I
      IY = J
      IG = 1
      GO TO 10
C PLOT POINT
  280 CALL DRAWG(' ',AX,AY,ICTYPE,LWTYPE,0,3)
      IX = I
      IY = J
      GO TO 10
C SET SCALES FOR FIRST PLOT
  290 CALL STARTG
      CALL GRPARM(3,LX,LY,MINX,MINY,LENX,LENY,ICH,ICW,ISLB,ISTCX,ISTCY,N
     1TCX,NTCY,IGSTYL)
      XS = FLOAT(LX - 1)/1023.
      YS = FLOAT(LY - 1)/780.
      ICX = 14.*XS + .5
      ICY = 22.*YS + .5
      IHY = 767.*YS + .5
      IBY = IZ
      IX = IBY
      IY = IHY
      LTYPE = LT(1)
      ICTYPE = ICOLOR(1)
      LWTYPE = LW(1)
      CALL INITGR(LX,LY,ICY,ICX)
      GO TO 10
C LAST PLOT
  300 NF = NF + 1
      IF (NF.GT.ND) GO TO 301
      CALL DRAWG(' ',0.,0.,0,0,0,0)
  301 IF (ID.EQ.1) CALL LABELF(NF)
      CALL DRAWG(' ',0.,0.,0,0,0,6)
      ND = ND + 1
      IF (ID.GT.1) GO TO 305
      CALL READN(IRC,NRC)
      IF (IRC.NE.1) GO TO 303
      CALL QUITG
      WRITE (6,93) NF
      STOP 1
  303 IF (IRC.EQ.2) ND = NRC - 1
      IF (NF.LE.ND) GO TO 305
      NF = 0
      IRD = 1
      IRB = 1
      IF (NF.LT.ND) GO TO 304
      CALL CHARSZ(ICY,ICX)
      CALL DRAWG(' ',0.,0.,0,0,0,0)
      LTYPE = LT(1)
      ICTYPE = ICOLOR(1)
      LWTYPE = LW(1)
  304 IG = 0
      LS = 0
      IF = 0
      IX = IZ
      IY = IHY
      GO TO 310
  305 CALL QUITG
      WRITE (6,93) NF
  310 RETURN
      END
*
*****************************************************
* LABELF -- PUTS LABEL IN GRAPHS FOR RASTER PLOTS
*****************************************************
*
      SUBROUTINE LABELF(NF)
C THIS SUBROUTINE PUTS LABEL IN GRAPHS FOR RASTER PLOTS
      CHARACTER*9 LBL
      SAVE LW,IC
      DATA LW,IC /1,7/
      CALL GRPARM(3,LX,LY,MINX,MINY,LENX,LENY,ICH,ICW,ISLB,ISTCX,ISTCY,N
     1TCX,NTCY,IGSTYL)
      CALL CHARSZ(ICH,ICW)
      AX = FLOAT(MINX + LENX - 9*ICW - 1)
      AY = FLOAT(MINY)
      IS = ICHAR('0')
      ID = 0
      N = NF
   10 ID = ID + 1
      N = N/10
      IF (N.GT.0) GO TO 10
      LBL = '#        '
      NT = NF
      LS = 10**(ID - 1)
      DO 20 I = 1, ID
      I1 = I + 1
      N = NT/LS
      LBL (I1:I1) = CHAR(N+IS)
      NT = NT - N*LS
      LS = LS/10
   20 CONTINUE
      N = ID + 1
      CALL DRSTRG(LBL,AX,AY,IC,LW,ICW,N)
      RETURN
      END
*
*****************************************************
* READN -- READS CHARACTERS
*****************************************************
*
      SUBROUTINE READN(IRC,NRC)
C THIS SUBROUTINE READS CHARACTERS AND OUTPUTS CODE AND VALUE
      CHARACTER*1 C
      CHARACTER*8 CX
      CHARACTER*40 CHR
   91 FORMAT (1X,A40)
      DATA CHR /' Q=QUIT, #=DISPLAY FRAME #, CR=CONTINUE '/
      IRC = 0
      NRC = 0
   10 CALL CINPUT(CX)
      C = CX(1:1)
      IF (C.NE.'?') GO TO 20
      CALL CLEAR
      WRITE (6,91) CHR
      GO TO 10
   20 IF ((C.EQ.'Q').OR.(C.EQ.'q')) IRC = 1
      IF (IRC.EQ.1) GO TO 30
      CALL EVALC(CX,IVAL,VAL,ID)
      IF (ID.EQ.0) GO TO 30
      NRC = IVAL
      IRC = 2
   30 RETURN
      END
*
*****************************************************
* SELFMP -- INITIALIZES DEFAULT MAPPING
*****************************************************
*
      SUBROUTINE SELFMP(IRX,IRY)
C THIS SUBROUTINE INITIALIZES DEFAULT MAPPING
      COMMON /RCOM/ MINX,LENX,MINY,LENY,ICX,ICY,
     1XMIN,XMAX,YMIN,YMAX,DX,DY,CSX,CSY
      MINX = 0
      LENX = IRX - 1
      MINY = 0
      LENY = IRY - 1
      XMIN = FLOAT(MINX)
      XMAX = FLOAT(MINX + LENX)
      YMIN = FLOAT(MINY)
      YMAX = FLOAT(MINY + LENY)
      DX = FLOAT(LENX)/(XMAX - XMIN)
      DY = FLOAT(LENY)/(YMAX - YMIN)
      RETURN
      END
*
*****************************************************
* MAPWIN -- SETS UP VARIABLE MAPPING
*****************************************************
*
      SUBROUTINE MAPWIN(XMN,XMX,YMN,YMX,MNX,LNX,MNY,LNY)
C THIS SUBROUTINE SETS UP VARIABLE MAPPING
      COMMON /RCOM/ MINX,LENX,MINY,LENY,ICX,ICY,
     1XMIN,XMAX,YMIN,YMAX,DX,DY,CSX,CSY
      MINX = MNX
      LENX = LNX
      MINY = MNY
      LENY = LNY
      XMIN = XMN
      XMAX = XMX
      YMIN = YMN
      YMAX = YMX
      DX = FLOAT(LENX)/(XMAX - XMIN)
      DY = FLOAT(LENY)/(YMAX - YMIN)
      RETURN
      END
*
*****************************************************
* CHARSZ -- SETS UP VARIABLE MAPPING
*****************************************************
*
      SUBROUTINE CHARSZ(ICH,ICW)
C THIS SUBROUTINE SETS UP VARIABLE MAPPING
      COMMON /RCOM/ MINX,LENX,MINY,LENY,ICX,ICY,
     1XMIN,XMAX,YMIN,YMAX,DX,DY,CSX,CSY
      CSX = FLOAT(ICW)/14.
      CSY = FLOAT(ICH)/16.
      RETURN
      END
*
*****************************************************
* ZIMAGE ZEROES OUT IMAGE
*****************************************************
*
      SUBROUTINE ZIMAGE(G,BLANK,LX,LY)
C THIS SUBROUTINE ZEROES OUT IMAGE
      CHARACTER*1 G(LX,LY)
      CHARACTER*1 BLANK
      DO 20 K = 1, LY
      DO 10 J = 1, LX
      G(J,K) = BLANK
   10 CONTINUE
   20 CONTINUE
      RETURN
      END
*
*****************************************************
* CDRAW -- DRAWS CHARACTER WITH ASCII CODE IC AT IX,IY
*****************************************************
*
      SUBROUTINE CDRAW(G,CTYPE,BLANK,LX,LY,IC,IX,IY,ICX,ICY,SX,SY,LWTYPE
     1)
C THIS SUBROUTINE DRAWS CHARACTER WITH ASCII CODE IC AT LOCATION IX,IY
C WITH SCALING FACTORS SX,SY
C FOR RASTER FILE
      CHARACTER*1 G(LX,LY)
      CHARACTER*1 CTYPE,BLANK
      DIMENSION LB(8)
      DIMENSION ICFLEN(94), ICFLOC(94)
      DIMENSION ICFON1(114), ICFON2(114), ICFONT(228)
      EQUIVALENCE (ICFON1(1), ICFONT(1)), (ICFON2(1), ICFONT(115))
      SAVE LW,NW,ICFLEN,ICFLOC,ICFONT
      DATA LW,NW /8,-2147483647/
      DATA ICFLEN /10,10,20,35,43,25,5,9,9,15,10,7,5,5,5,32,12,21,28,9,1
     19,23,9,35,23,10,12,7,10,7,18,47,13,23,25,16,14,9,32,12,15,16,9,7,8
     2,6,27,12,32,14,25,10,13,11,24,7,16,9,9,5,9,7,2,5,24,19,17,22,21,18
     3,28,13,10,11,9,11,24,13,19,22,22,11,21,18,16,11,24,7,22,9,23,10,23
     4,13/
      DATA ICFLOC /1,3,5,8,13,19,23,24,26,28,30,32,33,34,35,36,40,42,45,
     149,51,54,57,59,64,67,69,71,72,74,75,78,84,86,89,93,95,97,99,103,10
     25,107,109,111,112,113,114,118,120,124,126,130,132,134,136,139,140,
     3142,144,146,147,149,150,151,152,155,158,161,164,167,170,174,176,17
     48,180,182,184,187,189,192,195,198,200,203,206,208,210,213,214,217,
     5219,222,224,227/
      DATA ICFON1 /-184213676,1543503872,-225259652,2030043136,-21785408
     14,1895448655,145227776,-266266598,977822304,-2137868103,-135079911
     27,62914560,-266682549,1008470793,406341963,-164148918,974718726,37
     32244480,-100622159,-1010518880,1612842506,1073741824,-174735360,-1
     467562346,-1073741824,-200853612,-1073741824,-266682532,1358565552,
     5-261460132,1342177280,-217766608,-244752384,-184217600,-266686464,
     6-264207692,-959923574,1209402370,273613227,-224017137,545259520,-2
     757767222,-1433902496,1074397184,-257767222,-1433902225,-2036030848
     8,537001984,-133644182,1610612736,-88031128,1783244801,1048576,-896
     903392,-1608382454,709386336,-255146715,0,-121454432,-2107086262,67
     A1219744,1114139274,-1463812096,-94216064,-1563899222,671219744,-17
     B4747820,1392508928,-174747820,1395720192,-134059840,-259354716,671
     C08864,-234331456,-257767222,-1433055407,1342177280,-145258860,-182
     D0109754,1196968266,1519957700,-1028616142,335939600,123512736,-261
     E464064,210545320,-2046363542,1244135424,-90655036,-1028620222,3359
     F39610,805306368,209492648,-1533906944,212660327,1872756736,2126621
     G12,1610612736,-90655036,-1028620222,335939610,974103395,217082479,
     H-1398800384,-217641136,1559480256,-266205688,684682412,217759850,0
     I,217710592,207137952,211856384,-264207692/
      DATA ICFON2 /-959923574,1209402370,272629760,210545320,-2046427136
     1,-264207692,-959923574,1209402370,273638560,210545320,-2046386176,
     2-267319286,709386848,-2136815158,-1342177280,-184168692,-140928614
     34,-255839736,171622400,-255830774,1522532352,-255839741,86335314,1
     4887478700,-1393505792,-255814299,257337772,-1408435712,0,-18869452
     54,0,-256241664,-154482170,0,-265979360,-1610612736,-171483136,-244
     6152182,1779409956,35684513,217084552,-1972754430,1048576,-92700032
     7,1612843274,268435456,-99954777,-2010642942,545300736,-263566744,-
     82105515998,134811648,1089602227,-976830715,1157627904,-89603392,-1
     9604171702,1605018240,268500992,217084552,-1972764672,-184184997,15
     A43503872,-137943807,1048576,217743434,0,-238894556,83886080,149975
     B683,-2056974506,2022221472,149975688,-1972764672,-266313080,-19727
     C54430,2097152,-268382453,747416230,-2078014208,-99954773,-19432709
     D06,612672768,149963895,-1990197248,-250476534,675430498,-198741606
     E4,-205318906,135868168,1744830464,-91615327,-2145385976,-260025078
     F,1518338048,-260034045,86327122,1887478696,-1460631040,-87414783,2
     G034694,612672768,-259358710,0,-137968204,-1804377004,890636032,-17
     H1602092,1342177280,-205208138,-1770559914,890503936,-261979257,121
     I2833792/
      IS = IC - 32
      IF ((IS.LT.1).OR.(IS.GT.94)) GO TO 60
      ID = 0
      IF ((IS.EQ.71).OR.(IS.EQ.74).OR.(IS.EQ.80).OR.(IS.EQ.81).OR.(IS.EQ
     1.89)) ID = 4
      MW = NW - 1
      LWM = LW - 1
      LWP = LW + 1
      LEN = ICFLEN(IS)
      LENW = (LEN - 1)/LW + 1
      IOFF = ICFLOC(IS) - 1
      IL = LW
      MY = 0
      IM = 0
      CALL DMOVE(IX,IY,ICX,ICY)
C MAIN LOOP
      DO 50 I = 1, LENW
      IT = ICFONT(I+IOFF)
C     IF (IT.LT.0) IT = 4294967296 + IT
      LB(1) = IT
      IF (LB(1).LT.0) LB(1) = LB(1) - MW
C DECODE COORDINATES
      DO 10 J = 1, LWM
      LB(J+1) = LB(J)/16
      LB(J) = LB(J) - LB(J+1)*16
   10 CONTINUE
      IF (IT.LT.0) LB(LW) = LB(LW) + 8
      IF (I.EQ.LENW) IL = LEN - (I - 1)*LW
C DRAW CHARACTER
      DO 40 J = 1, IL
      IT = LB(LWP-J)
      IF (MY.EQ.1) GO TO 30
C X COORDINATE OR MOVE FLAG
      IF (IT.EQ.15) GO TO 20
      IT = SX*FLOAT(IT) + .5
      JX = IX + IT
      MY = 1
      GO TO 40
   20 IM = 1
      GO TO 40
C Y COORDINATE
   30 IT = SY*FLOAT(IT - ID) + .5
      JY = IY + IT
      IF (IM.EQ.0) CALL DLINE(G,CTYPE,BLANK,LX,LY,JX,JY,ICX,ICY,LWTYPE)
      IF (IM.EQ.1) CALL DMOVE(JX,JY,ICX,ICY)
      MY = 0
      IM = 0
   40 CONTINUE
   50 CONTINUE
      CALL DMOVE(IX,IY,ICX,ICY)
   60 RETURN
      END
*
*****************************************************
* DMOVE -- MOVES CURSOR
*****************************************************
*
      SUBROUTINE DMOVE(I,J,ICX,ICY)
C THIS SUBROUTINE MOVES CURSOR TO (I,J)
      ICX = I
      ICY = J
      RETURN
      END
*
*****************************************************
* DLINE -- DRAWS LINE
*****************************************************
*
      SUBROUTINE DLINE(G,CTYPE,BLANK,LX,LY,I,J,ICX,ICY,LWTYPE)
C THIS SUBROUTINE DRAWS LINE FROM (ICX,ICY) TO (I,J)
C WITH CTYPE CHARACTER
      CHARACTER*1 G(LX,LY)
      CHARACTER*1 CTYPE,BLANK
      CHARACTER*1 XX
      XX = CTYPE
      IF (LWTYPE.LT.1) GO TO 60
      II = I - ICX
      JJ = J - ICY
      KK = 1
      AJI = 1.
      ALWTYPE = .5*FLOAT(LWTYPE)
      IF (IABS(JJ).GT.IABS(II)) GO TO 30
      IF (II.LT.0) KK = -1
      IF (II.NE.0) AJI = FLOAT(JJ)/FLOAT(II)
      DO 20 N = 1, LWTYPE
      AJ = FLOAT(ICY + N) - ALWTYPE
      DO 10 L = ICX, I, KK
      IF ((L.LT.0).OR.(L.GE.LX)) GO TO 10
      M = AJI*FLOAT(L - ICX) + AJ
      IF ((M.LT.0).OR.(M.GE.LY)) GO TO 10
      IF (G(L+1,M+1).EQ.BLANK) G(L+1,M+1) = CTYPE
      IF (G(L+1,M+1).NE.CTYPE) G(L+1,M+1) = XX
C     IF (G(L+1,M+1).NE.CTYPE) G(L+1,M+1) = CHAR(IOR(ICHAR(G(L+1,M+1)),I
C    1CHAR(CTYPE)))
   10 CONTINUE
   20 CONTINUE
      GO TO 60
   30 IF (JJ.LT.0) KK = -1
      AJI = FLOAT(II)/FLOAT(JJ)
      DO 50 N = 1, LWTYPE
      AI = FLOAT(ICX + N) - ALWTYPE
      DO 40 M = ICY, J, KK
      IF ((M.LT.0).OR.(M.GE.LY)) GO TO 40
      L = AJI*FLOAT(M - ICY) + AI
      IF ((L.LT.0).OR.(L.GE.LX)) GO TO 40
      IF (G(L+1,M+1).EQ.BLANK) G(L+1,M+1) = CTYPE
      IF (G(L+1,M+1).NE.CTYPE) G(L+1,M+1) = XX
C     IF (G(L+1,M+1).NE.CTYPE) G(L+1,M+1) = CHAR(IOR(ICHAR(G(L+1,M+1)),I
C    1CHAR(CTYPE)))
   40 CONTINUE
   50 CONTINUE
   60 ICX = I
      ICY = J
      RETURN
      END
*
*****************************************************
* DPNT -- DRAWS POINT WITH CTYPE CHARACTER INTO G
*****************************************************
*
      SUBROUTINE DPNT(G,CTYPE,BLANK,LX,LY,I,J,ICX,ICY,LWTYPE)
C THIS SUBROUTINE DRAWS POINT AT (I,J) WITH CTYPE CHARACTER
      CHARACTER*1 G(LX,LY)
      CHARACTER*1 CTYPE,BLANK
      CHARACTER*1 XX
      XX = CTYPE
      IF (LWTYPE.LT.1) GO TO 30
      IS = (LWTYPE + 1)/2
      N1 = I - IS
      K1 = J - IS
      DO 20 K = 1, LWTYPE
      M = K1 + K
      IF ((M.LT.0).OR.(M.GE.LY)) GO TO 20
      DO 10 N = 1, LWTYPE
      L = N1 + N
      IF ((L.LT.0).OR.(L.GE.LX)) GO TO 10
      IF (G(L+1,M+1).EQ.BLANK) G(L+1,M+1) = CTYPE
      IF (G(L+1,M+1).NE.CTYPE) G(L+1,M+1) = XX
C     IF (G(L+1,M+1).NE.CTYPE) G(L+1,M+1) = CHAR(IOR(ICHAR(G(L+1,M+1)),I
C    1CHAR(CTYPE)))
   10 CONTINUE
   20 CONTINUE
   30 ICX = I
      ICY = J
      RETURN
      END
*
*****************************************************
* DASHLN -- DRAWS A DASHED LINE
*****************************************************
*
      SUBROUTINE DASHLN(G,CTYPE,BLANK,LX,LY,I,J,ICX,ICY,L,LWTYPE)
C THIS SUBROUTINE DRAWS A DASHED LINE FROM (ICX,ICY) TO (I,J)
      CHARACTER*1 G(LX,LY)
      CHARACTER*1 CTYPE,BLANK
      DIMENSION LT(4,4)
      SAVE NS,NL,LTYPE
      DATA LTYPE /0/
      DATA LT /5,5,5,5,14,6,4,6,9,6,9,6,23,7,23,7/
      IF ((L.EQ.LTYPE).OR.(L.LT.0).OR.(L.GT.7)) GO TO 10
      LTYPE = L
      IF ((LTYPE.EQ.0).OR.(LTYPE.GT.4)) GO TO 50
      NS = 0
      NL = LT(NS+1,LTYPE)
      GO TO 20
   10 IF ((LTYPE.EQ.0).OR.(LTYPE.GT.4)) GO TO 50
C SOFTWARE DASHED LINE
   20 IX0 = ICX
      IY0 = ICY
      COST = FLOAT(I - ICX)
      SINT = FLOAT(J - ICY)
      ALEN = SQRT(COST*COST + SINT*SINT)
      LEN = ALEN + .5
      IF (NL.GE.LEN) GO TO 40
      COST = COST/ALEN
      SINT = SINT/ALEN
   30 ANL = FLOAT(NL)
      IX = FLOAT(ICX) + ANL*COST + .5
      IY = FLOAT(ICY) + ANL*SINT + .5
      IT = NS - (NS/2)*2
      IF (IT.EQ.0) CALL DLINE(G,CTYPE,BLANK,LX,LY,IX,IY,IX0,IY0,LWTYPE)
      IF (IT.EQ.1) CALL DMOVE(IX,IY,IX0,IY0)
      NS = NS + 1
      IF (NS.EQ.4) NS = 0
      NL = NL + LT(NS+1,LTYPE)
      IF (NL.LT.LEN) GO TO 30
   40 IT = NS - (NS/2)*2
      IF (IT.EQ.0) CALL DLINE(G,CTYPE,BLANK,LX,LY,I,J,IX0,IY0,LWTYPE)
      IF (IT.EQ.1) CALL DMOVE(I,J,IX0,IY0)
      NL = NL - LEN
      IF (NL.GT.0) GO TO 60
      NS = NS + 1
      IF (NS.EQ.4) NS = 0
      NL = LT(NS+1,LTYPE)
      GO TO 60
C SOLID LINE
   50 CALL DLINE(G,CTYPE,BLANK,LX,LY,I,J,ICX,ICY,LWTYPE)
   60 ICX = I
      ICY = J
      RETURN
      END
*
*****************************************************
* WIMAGE -- WRITES IMAGE (G) TO TERMINAL
*****************************************************
*
      SUBROUTINE WIMAGE(G,LX,LY)
C THIS SUBROUTINE WRITES IMAGE TO TERMINAL
      CHARACTER*1 G(LX,LY)
   91 FORMAT (1H1,128A1)
   92 FORMAT (1X,128A1)
      LY1 = LY + 1
      DO 10 K = 1, LY
      K1 = LY1 - K
      IF (K.EQ.1) WRITE (6,91) (G(J,K1),J=1,LX)
      IF (K.GT.1) WRITE (6,92) (G(J,K1),J=1,LX)
   10 CONTINUE
      RETURN
      END
*
*****************************************************
* CWRITE -- WRITES A CHARACTER INTO G
*****************************************************
*
      SUBROUTINE CWRITE(G,C,LX,LY,IX,IY,ICX,ICY)
C THIS SUBROUTINE WRITES CHARACTER AT LOCATION (IX,IY)
      CHARACTER*1 G(LX,LY)
      CHARACTER*1 C
      G(IX+1,IY+1) = C
      ICX = IX
      ICY = IY
      RETURN
      END
*
*****************************************************
* header -- writes header in movie file
*****************************************************
*
      SUBROUTINE HEADER(IFRMT,LX,LY,NBIT)

C THIS SUBROUTINE WRITES HEADER IN MOVIE FILE
C IF IFRMT = 1,2,3  UCLA FORMAT HEADER = IFRMT, LX, LY, NBIT
C LX, LY IS THE SIZE OF THE IMAGE, AND NBIT IS THE NUMBER OF COLORS
C IF IFRMT = 4,   MFE FORMAT HEADER = FRMT
C WHERE FRMT IS A 1 BYTE CHARACTER VARIABLE GIVEN BY:
C 241 FOR CGA, 242 FOR EGA, 243 FOR VGA
C IF IFRMT = 0, NO HEADER

      CHARACTER*1 CHR(16)
      DIMENSION IMGSIZ(4)
      IF ((IFRMT.LT.1).OR.(IFRMT.GT.4)) GO TO 20
      IF (IFRMT.EQ.4) GO TO 10
      IMGSIZ(1) = IFRMT
      IMGSIZ(2) = LX
      IMGSIZ(3) = LY
      IMGSIZ(4) = NBIT
C UCLA FORMAT HEADER
      CALL CONVIC(IMGSIZ,CHR,16,4)
      CALL BUFFWR(CHR,16,IFRMT)
      GO TO 20
C MFE FORMAT HEADER
   10 CHR(1) = CHAR(0)
      IF ((LX.EQ.320).AND.(LY.EQ.200).AND.(NBIT.EQ.2)) CHR(1)=CHAR(241)
      IF ((LX.EQ.640).AND.(LY.EQ.350).AND.(NBIT.EQ.1)) CHR(1)=CHAR(242)
      IF ((LX.EQ.320).AND.(LY.EQ.200).AND.(NBIT.EQ.8)) CHR(1)=CHAR(243)
      CALL BUFFWR(CHR,1,IFRMT)
   20 RETURN
      END
*
*****************************************************
* wrpal -- write palette for mve vga format files
*****************************************************
*
      SUBROUTINE WRPAL(PAL,NPAL,IFRMT)
C THIS SUBROUTINE WRITES PALETTE FOR MFE VGA FORMAT FILES
      CHARACTER*1 PAL(768)
      CHARACTER*1 COLOR(24)
      CHARACTER*1 CHR
      DIMENSION ICOLOR(24)
      SAVE ICOLOR
      DATA ICOLOR /0,0,0,0,0,1,0,1,0,0,1,1,1,0,0,1,0,1,1,1,0,1,1,1/
      IF (IFRMT.NE.4) GO TO 80
C WRITE DEFAULT PALETTE FOR VGA MODE
      CHR = CHAR(240)
      CALL BUFFWR(CHR,1,IFRMT)
      IF (NPAL.GT.0) GO TO 40
         DO 10 I = 1, 24
            COLOR(I) = CHAR(63*ICOLOR(I))
   10    CONTINUE
      CALL BUFFWR(COLOR,24,IFRMT)
         DO 20 I = 1, 24
            COLOR(I) = CHAR(63)
   20    CONTINUE
         DO 30 I = 1, 31
            CALL BUFFWR(COLOR,24,IFRMT)
   30    CONTINUE
      GO TO 80
C WRITE USER PALETTE FOR VGA MODE
   40 LEN = 3*NPAL
      IF (LEN.GT.768) LEN = 768
      CALL BUFFWR(PAL,LEN,IFRMT)
      IF (LEN.EQ.768) GO TO 80
C PAD PALETTE WITH WHITE
      NL = (768 - LEN)
      N = NL/24
      NL = NL - 24*N
      DO 50 I = 1, 24
      COLOR(I) = CHAR(63)
   50 CONTINUE
C     write( 6,* ) ' n = ', n
      IF (N.EQ.0) GO TO 70
      DO 60 I = 1, N
      CALL BUFFWR(COLOR,24,IFRMT)
   60 CONTINUE
   70 IF (NL.GT.0) CALL BUFFWR(COLOR,NL,IFRMT)
   80 RETURN
      END
*
*****************************************************
* GIMAGE -- COMPRESSES CHARACTER IMAGE DATA INTO CHARACTER DATA
*****************************************************
*
      SUBROUTINE GIMAGE (G,IMAGE,LX,LY,LZ,NBIT,INTRL)
C THIS SUBROUTINE COMPRESSES CHARACTER IMAGE DATA INTO CHARACTER DATA
C NBIT = NUMBER OF BITS PER PIXEL
C LZ = (LX - 1)/NPIX + 1, WHERE NPIX = 8/NBIT
      CHARACTER*1 G(LX,LY)
      CHARACTER*1 BLANK
      CHARACTER*1 IMAGE(LZ,LY)
C     DATA BLANK /' '/
      BLANK = CHAR(0)
      NPIX = 8/NBIT
      NTC = 2**NBIT
      LY1 = LY + 1
      LYH = (LY - 1)/2 + 1
      DO 40 K = 1, LY
         K1 = LY1 - K
         K2 = K
         K3 = (K - 1)/2
         IF (INTRL.EQ.1) K2 = LYH*((K2 - 1) - 2*K3) + K3 + 1
         DO 30 J = 1, LZ
            J1 = (J - 1)*NPIX
            ITC = 0
            DO 20 I = 1, NPIX
               IT = 0
               J2 = J1 + I
               IF (J2.GT.LX) GO TO 10
C     IF (G(J2,K1).NE.BLANK) IT = 1
               IS = ICHAR(G(J2,K1))
               IT = IS - (IS/NTC)*NTC
   10          ITC = NTC*ITC + IT
   20       CONTINUE
            IMAGE(J,K2) = CHAR(ITC)
   30    CONTINUE
   40 CONTINUE
      RETURN
      END
*
*****************************************************
* GRIMAGE -- COMPRESSES CHARACTER IMAGE DATA INTO CHARACTER DATA
* THIS ROUTINE COMPRESSES THE STUFF FROM G INTO IMAGE
*****************************************************
*
      SUBROUTINE GRIMAGE (G,IMAGE,NX,NY,NXV,LX,LY,LZ,NBIT,INTRL,IREV)
C THIS SUBROUTINE COMPRESSES CHARACTER IMAGE DATA INTO CHARACTER DATA
C IF LX < NX OR LY < NY, THEN IMAGE IS TRUNCATED
C IF LX > NX OR LY > NY, THEN IMAGE IS PADDED WITH NULLS
C NBIT = NUMBER OF BITS PER PIXEL
C LZ = (LX - 1)/NPIX + 1, WHERE NPIX = 8/NBIT
C INTRL = (0,1) = (NO,YES) INTERLACE IMAGE
C IREV = (0,1) = (NO,YES) FLIP IMAGE VERTICALLY
      CHARACTER*1 G(NXV,NY)
      CHARACTER*1 BLANK
      CHARACTER*1 IMAGE(LZ,LY)
C     DATA BLANK /' '/
      BLANK = CHAR(0)
      NPIX = 8/NBIT
      NTC = 2**NBIT
      LY1 = LY + 1
      LYH = (LY - 1)/2 + 1
      IF (NY.LT.LY) LY1 = NY + 1
C LOOP OVER ROWS
      DO 60 K = 1, LY
      K1 = K
      IF (IREV.EQ.1) K1 = LY1 - K
      K2 = K
      K3 = (K - 1)/2
      IF (INTRL.EQ.1) K2 = LYH*((K2 - 1) - 2*K3) + K3 + 1
      IF ((K1.LT.1).OR.(K1.GT.NY)) GO TO 40
C LOOP OVER COLUMNS
      DO 30 J = 1, LZ
      J1 = (J - 1)*NPIX
      ITC = 0
C EXTRACT LOW ORDER NBITS FROM G ARRAY
      DO 20 I = 1, NPIX
      IT = 0
      J2 = J1 + I
      IF ((J2.GT.NX).OR.(J2.GT.LX)) GO TO 10
C     IF (G(J2,K1).NE.BLANK) IT = 1
      IS = ICHAR(G(J2,K1))
      IT = IS - (IS/NTC)*NTC
   10 ITC = NTC*ITC + IT
   20 CONTINUE
      IMAGE(J,K2) = CHAR(ITC)
   30 CONTINUE
      GO TO 60
C PAD Y VALUES WITH NULLS
   40 ITC = 0
      DO 50 J = 1, LZ
      IMAGE(J,K2) = CHAR(ITC)
   50 CONTINUE
   60 CONTINUE
      RETURN
      END
*
*****************************************************
* PCTSAV -- COMPRESSES SUCCESSIVE IMAGES
*****************************************************
*
      SUBROUTINE PCTSAV(IMAGE,JMAGE,IMG,LINE,LIMG,IXOR,LZ,LY,LENB,LENG,L
     1ZG,IFRMT,INTRL)
C THIS SUBROUTINE PERFORMS COMPRESSION OF SUCCESSIVE IMAGES.
C INPUT IS IN ARRAY IMAGE, AND COMPRESSED OUTPUT IS IN ARRAY IMG.
C IF IXOR = 1, THEN THE CURRENT IMAGE IS XORED WITH THE PREVIOUS IMAGE
C BEFORE COMPRESSION, AND THE CURRENT IMAGE IS SAVE IN THE ARRAY JMAGE.
C LINE AND LIMG ARE SCRATCH ARRAYS NEEDED BY SUBROUTINE COMPRS
      CHARACTER*1 IMAGE(LENB), JMAGE(LENB)
      CHARACTER*1 IMG(LENG), LINE(LZ), LIMG(LZG)
      CHARACTER*1 CHR(4)
   91 FORMAT(36H COMPRESSED IMAGE OVERFLOW, IBMAX = ,I20,8H LENG = ,I20)
C COMPRESS IMAGE
      IF (INTRL.EQ.1) GO TO 10
      IF (IXOR.EQ.1) GO TO 20
      IF (IFRMT.EQ.3) CALL COPYIM(IMAGE,IMG,LENB,LENG,IBMAX)
      IF (IFRMT.NE.3) CALL COMPRS(IMAGE,IMG,LINE,LIMG,LZ,LY,LZ,LENB,LENG
     1,LZG,IBMAX)
      IF (IBMAX.GT.LENG) WRITE (6,91) IBMAX, LENG
      GO TO 50
C INTERLACED IMAGE
   10 LYH = LY/2
      LYH1 = (LY - 1)/2 + 1
      LEN = LZ*LYH1
      CALL COMPRS(IMAGE,IMG,LINE,LIMG,LZ,LYH1,LZ,LEN,LENG,LZG,IBMAX)
      IMG(IBMAX) = CHAR(240)
      JBMAX = IBMAX
      IT1 = LEN + 1
      IT2 = IBMAX + 1
      LEN = LZ*LYH
      CALL COMPRS(IMAGE(IT1),IMG(IT2),LINE,LIMG,LZ,LYH,LZ,LEN,LENG,LZG,I
     1BMAX)
      IBMAX = IBMAX + JBMAX
      IF (IBMAX.GT.LENG) WRITE (6,91) IBMAX, LENG
      GO TO 60
C XOR AND COMPRESS SUCCESSIVE IMAGES
   20 DO 30 I = 1, LENB
      JMAGE(I) = CHAR(IEOR(ICHAR(IMAGE(I)),ICHAR(JMAGE(I))))
C     JMAGE(I) = CHAR(ICHAR(IMAGE(I)).XOR.ICHAR(JMAGE(I)))
   30 CONTINUE
      IF (IFRMT.EQ.3) CALL COPYIM(JMAGE,IMG,LENB,LENG,IBMAX)
      IF (IFRMT.NE.3) CALL COMPRS(JMAGE,IMG,LINE,LIMG,LZ,LY,LZ,LENB,LENG
     1,LZG,IBMAX)
      IF (IBMAX.GT.LENG) WRITE (6,91) IBMAX, LENG
C SAVE OLD IMAGE
      DO 40 I = 1, LENB
      JMAGE(I) = IMAGE(I)
   40 CONTINUE
C WRITE COMPRESSED IMAGE TO DISK
   50 IF (IFRMT.EQ.4) GO TO 60
      CALL CONVIC(IBMAX,CHR,4,1)
      CALL BUFFWR(CHR,4,IFRMT)
   60 CALL BUFFWR(IMG,IBMAX,IFRMT)
      RETURN
      END
*
*****************************************************
* COPYIM -- MOVES CHARACTER DATA FROM IMAGE TO IMG
*****************************************************
*
      SUBROUTINE COPYIM (IMAGE,IMG,LMAX,LMAXG,IB)
C THIS SUBROUTINE MOVES CHARACTER DATA FROM ARRAY IMAGE TO IMG
      CHARACTER*1 IMAGE(LMAX), IMG(LMAXG)
      IB = LMAX
         DO 10 I = 1, IB
            IMG(I) = IMAGE(I)
   10    CONTINUE
      RETURN
      END
*
*****************************************************
* COMPRS -- COMPRESSES BINARY DATA
*****************************************************
*
      SUBROUTINE COMPRS (IMAGE,IMG,LINE,LIMG,LENL,NUML,LENV,LMAX,LMAXG,L
     1ENLG,IB)
C THIS SUBROUTINE COMPRESSES BINARY DATA, USING THE ALGORITHM
C BY R. H. FROBOSE, JR., LAWRENCE LIVERMORE LAB REPORT UCRL-51858
C WRITTEN BY VIKTOR K. DECYK, UCLA
C INPUT IS IN ARRAY IMAGE, OF SIZE LMAX = LENV*NUML, WHERE
C LENV = SPACING BETWEEN ROWS IN BYTES, NUML = NUMBER OF LINES
C LENL = LENGTH OF LINE IN BYTES
C OUTPUT IS IN ARRAY IMG, OF MAXIMUM SIZE LMAXG = LMAX+(LMAX-1)/3+1,
C AND ACTUAL SIZE GIVEN IN VARIABLE IB
C LINE AND LIMG ARE SCRATCH ARRAYS, OF SIZE LENL AND LENLG RESPECTIVELY,
C WHERE LENLG = LENL+(LENL-1)/3+1
      CHARACTER*1 IMAGE(LMAX), IMG(LMAXG), LINE(LENL), LIMG(LENLG)
      INTEGER IC0,IC1,IC2,IC3
      SAVE ICN,IC0,IC1,IC2,IC3,ISTYLE
      DATA ICN,IC0,IC1,IC2,IC3 /64,0,64,128,192/
C ISTYLE = (0,1,2) = (NO XOR,XOR,BOTH) SUCCESSIVE LINES
      DATA ISTYLE /2/
      ICNR = ICN/2 - 2
C PREVENT XOR ON FIRST LINE
      IB = LMAX
      IB0 = 0
      IL = 0
      GO TO 120
   10 IF (IXOR.EQ.1) GO TO 100
      IF ((ISTYLE.GT.0).AND.(LIB.GT.(IB-IB0))) GO TO 30
C XORED LINE IS LONGER
      DO 20 J = 1, LIB
      IMG(J+IB0) = LIMG(J)
   20 CONTINUE
      IB = IB0 + LIB
C START NEW LINE
   30 IL = IL + LENV
      IF (IL.GE.LMAX) GO TO 450
C LOOKING FOR IDENTICAL LINES
      IRL = 0
      IL1 = IL - LENV
      I = 0
   40 IF (I.EQ.LENL) GO TO 50
      I = I + 1
      IF (IMAGE(I+IL).EQ.IMAGE(I+IL1)) GO TO 40
C LINE NOT IDENTICAL
      IF (IRL.EQ.0) GO TO 60
C REPEAT PREVIOUS LINE IRL TIMES.  DIFFERENT LINE FOLLOWS.
      IB = IB + 1
      ITC = IC3 + IRL
      IMG(IB) = CHAR(ITC)
      GO TO 60
C LINE IDENTICAL
   50 IRL = IRL + 1
      IL = IL + LENV
      I = 0
      IF ((IL.LT.LMAX).AND.(IRL.LT.ICNR)) GO TO 40
C REPEAT PREVIOUS LINE IRL TIMES.  BUFFER FULL.
      IB = IB + 1
      ITC = IC3 + IRL
      IMG(IB) = CHAR(ITC)
      IF (IL.GE.LMAX) GO TO 450
      IRL = 0
      I = 0
      GO TO 40
C TEST WHETHER TO SKIP XOR
   60 IB0 = IB
      IF (ISTYLE.EQ.0) GO TO 120
C  XOR CURRENT LINE
      IL1 = IL - LENV
      DO 70 I = 1, LENL
      LINE(I) = CHAR(IEOR(ICHAR(IMAGE(I+IL)),ICHAR(IMAGE(I+IL1))))
C     LINE(I) = CHAR(ICHAR(IMAGE(I+IL)).XOR.ICHAR(IMAGE(I+IL1)))
   70 CONTINUE
      LIB = 1
      ITC = IC3 + ICNR + 1
      LIMG(LIB) = CHAR(ITC)
      IXOR = 1
      GO TO 140
C SECOND PASS
C SAVE PREVIOUS COMPRESSED LINE
  100 DO 110 J = 1, LIB
      IMG(J+IB) = LIMG(J)
  110 CONTINUE
      IB = IB + LIB
C TEST WHETHER TO PERFORM XOR
      IF (ISTYLE.EQ.1) GO TO 30
  120 DO 130 I = 1, LENL
      LINE(I) = IMAGE(I+IL)
  130 CONTINUE
      LIB = 0
      IXOR = 0
  140 I = 1
      IZ = 0
      ID = 0
  150 IF (LINE(I).EQ.CHAR(0)) GO TO 400
C NON-ZERO BYTES
C LOOKING FOR DIFFERENT BYTES
  200 ID = ID + 1
      IF (I.EQ.LENL) GO TO 260
      IF (ID.EQ.ICN) GO TO 240
      I = I + 1
      IF (LINE(I).NE.CHAR(0)) GO TO 205
      IF ((I.EQ.LENL).OR.(LINE(I+1).NE.CHAR(0))) GO TO 200
      GO TO 220
  205 IF (LINE(I).NE.LINE(I-1)) GO TO 200
      IF (ID.EQ.1) GO TO 300
      IF ((I.EQ.LENL).OR.(LINE(I+1).NE.LINE(I))) GO TO 200
      ID = ID - 1
C PLOT NEXT ID BYTES AS THEY APPEAR. PAIR FOLLOWS.
      LIB = LIB + 1
      ITC = IC2 + ID - 1
      LIMG(LIB) = CHAR(ITC)
      I1 = (I - ID - 2)
      DO 210 J = 1, ID
      LIMG(J+LIB) = LINE(J+I1)
  210 CONTINUE
      LIB = LIB + ID
      GO TO 300
C PLOT NEXT ID BYTES AS THEY APPEAR. ZERO FOLLOWS.
  220 LIB = LIB + 1
      ITC = IC2 + ID - 1
      LIMG(LIB) = CHAR(ITC)
      I1 = (I - ID - 1)
      DO 230 J = 1, ID
      LIMG(J+LIB) = LINE(J+I1)
  230 CONTINUE
      LIB = LIB + ID
      ID = 0
      GO TO 400
C PLOT NEXT ID BYTES AS THEY APPEAR.  BUFFER FULL.
  240 LIB = LIB + 1
      ITC = IC2 + ID - 1
      LIMG(LIB) = CHAR(ITC)
      I1 = (I - ID)
      DO 250 J = 1, ID
      LIMG(J+LIB) = LINE(J+I1)
  250 CONTINUE
      LIB = LIB + ID
      I = I + 1
      ID = 0
      GO TO 150
C PLOT NEXT ID BYTES AS THEY APPEAR.  LINE FULL.
  260 LIB = LIB + 1
      ITC = IC2 + ID - 1
      LIMG(LIB) = CHAR(ITC)
      I1 = (I - ID)
      DO 270 J = 1, ID
      LIMG(J+LIB) = LINE(J+I1)
  270 CONTINUE
      LIB = LIB + ID
      GO TO 10
C NON-ZERO BYTES
C LOOKING FOR IDENTICAL BYTES
  300 IR = 1
  310 IR = IR + 1
      IF (I.EQ.LENL) GO TO 330
      IF (IR.EQ.ICN) GO TO 320
      I = I + 1
      IF (LINE(I).EQ.LINE(I-1)) GO TO 310
C REPEAT NEXT BYTE IR TIMES.  DIFFERENT BYTE FOLLOWS.
      LIB = LIB + 1
      ITC = IC0 + IR - 1
      LIMG(LIB) = CHAR(ITC)
      LIB = LIB + 1
      LIMG(LIB) = LINE(I-1)
      ID = 0
      GO TO 150
C REPEAT NEXT BYTE IR TIMES.  BUFFER FULL.
  320 LIB = LIB + 1
      ITC = IC0 + IR - 1
      LIMG(LIB) = CHAR(ITC)
      LIB = LIB + 1
      LIMG(LIB) = LINE(I)
      I = I + 1
      ID = 0
      GO TO 150
C REPEAT NEXT BYTE IR TIMES.  LINE FULL.
  330 LIB = LIB + 1
      ITC = IC0 + IR - 1
      LIMG(LIB) = CHAR(ITC)
      LIB = LIB + 1
      LIMG(LIB) = LINE(I)
      GO TO 10
C ZERO BYTES
  400 IZ = IZ + 1
      IF (I.EQ.LENL) GO TO 430
      I = I + 1
      IF (LINE(I).EQ.CHAR(0)) GO TO 400
  410 IF (IZ.LE.ICN) GO TO 420
C SKIP NEXT IZ BYTES, AND PLOT THE FOLLOWING BYTE.  BUFFER FULL.
      LIB = LIB + 1
      ITC = IC1 + ICN - 1
      LIMG(LIB) = CHAR(ITC)
      LIB = LIB + 1
      LIMG(LIB) = CHAR(0)
      IZ = IZ - (ICN + 1)
      IF (IZ.EQ.0) GO TO 200
      GO TO 410
C SKIP NEXT IZ BYTES, AND PLOT THE FOLLOWING BYTE.  NEXT BYTE NON-ZERO.
  420 LIB = LIB + 1
      ITC = IC1 + IZ - 1
      LIMG(LIB) = CHAR(ITC)
      LIB = LIB + 1
      LIMG(LIB) = LINE(I)
      IF (I.EQ.LENL) GO TO 10
      I = I + 1
      IZ = 0
      GO TO 150
C PLOT END-OF-LINE SENTINEL
  430 LIB = LIB + 1
      LIMG(LIB) = CHAR(IC3)
      GO TO 10
C END-OF-FRAME BYTE
  450 IB = IB + 1
      IMG(IB) = CHAR(0)
      RETURN
      END
*
*****************************************************
* BUFFWR -- PACKS IMAGE INTO BUFFER AND WRITES WHEN FULL
*****************************************************
*
      SUBROUTINE BUFFWR(LINE,N,IFRMT)
C THIS SUBROUTINE PACKS IMAGE DATA INTO BUFFER AND WRITES WHEN FULL
C INPUT IS IN CHARACTER ARRAY LINE, AND N CHARACTERS ARE TO BE WRITTEN
      CHARACTER*1 C0
      CHARACTER*1 LINE(*)
      CHARACTER*1 LOUT(80)
      DIMENSION IOUT(20)
      SAVE LEN,LMAX,LW,LOUT
      DATA LEN,LMAX,LW /0,80,4/
      L = (LMAX - 1)/LW + 1
      NC = N
      NCR = N
      I = 0
      IF (N.GT.0) GO TO 10
      IF (LEN.GE.0) GO TO 60
      GO TO 70
   10 NCR = (NC + LEN) - LMAX
      IF (NCR.GE.0) NC = LMAX - LEN
   20 IF (LEN.GT.0) GO TO 40
C     C0 = CHAR(0)
      C0 = CHAR(232)
      DO 30 J = 1, LMAX
      LOUT(J) = C0
   30 CONTINUE
   40 DO 50 J = 1, NC
      LOUT(J+LEN) = LINE(I+J)
   50 CONTINUE
      I = I + NC
      LEN = LEN + NC
      IF (NCR.LT.0) GO TO 70
   60 ITC = LEN
      IF (IFRMT.EQ.4) ITC = LMAX
      IF (LEN.GT.0) CALL PCTOUT(LOUT,IOUT,ITC,L,IFRMT)
      LEN = 0
      IF (N.EQ.0) CALL PCTOUT(LOUT,IOUT,LEN,L,IFRMT)
      IF (NCR.EQ.0) GO TO 70
      NC = NCR
      NCR = NCR - LMAX
      IF (NCR.GE.0) NC = LMAX
      GO TO 20
   70 RETURN
      END
*
*****************************************************
* PCTOUT -- WRITES COMPRESSED RASTER FILE IN APPROP FORM
*****************************************************
*
      SUBROUTINE PCTOUT(LOUT,IOUT,LEN,N,IFRMT)
C THIS SUBROUTINE WRITES COMPRESSED RASTER FILE IN APPROPRIATE FORM
      CHARACTER*1 LOUT(*)
      DIMENSION IOUT(N)
   91 FORMAT (80A1)
C IFRMT = 2 = ENCRYPT THE FILE
      IF (IFRMT.EQ.2) GO TO 10
      IF (LEN.GT.0) WRITE (19,91) (LOUT(J),J=1,LEN)
      GO TO 20
   10 CALL CONVCI(LOUT,IOUT,LEN,N)
      CALL ENCODE(IOUT,LEN)
   20 RETURN
      END
*
*****************************************************
* CONVIC -- CONVERTS PACKED INTEGER TO CHARACTER DATA
*****************************************************
*
      SUBROUTINE CONVIC(LIN,CHR,LEN,N)
C THIS SUBROUTINE CONVERTS PACKED INTEGER TO CHARACTER DATA
C SHOULD HAVE N = (LEN - 1)/LW + 1
C DIMENSION LB(LW)
C MW = -2**(8*LW-1)
      CHARACTER*1 CHR(LEN)
      DIMENSION LIN(N)
      DIMENSION LB(4)
      SAVE LW,NW
      DATA LW,NW /4,-2147483647/
      IF (LEN.LT.1) GO TO 70
      MW = NW - 1
      L = (LEN - 1)/LW + 1
      M = LEN/LW
      MR = LEN - M*LW
      LWM = LW - 1
      LWP = LW + 1
      IF (M.EQ.0) GO TO 40
      DO 30 I = 1, M
      I1 = (I - 1)*LW
      LB(1) = LIN(I)
      IF (LB(1).LT.0) LB(1) = LB(1) - MW
      DO 10 J = 1, LWM
      LB(J+1) = LB(J)/256
      LB(J) = LB(J) - LB(J+1)*256
   10 CONTINUE
      IF (LIN(I).LT.0) LB(LW) = LB(LW) + 128
      DO 20 J = 1, LW
      CHR(I1 + J) = CHAR(LB(LWP - J))
   20 CONTINUE
   30 CONTINUE
   40 IF (MR.EQ.0) GO TO 70
      I1 = M*LW
      LB(1) = LIN(L)
      IF (LB(1).LT.0) LB(1) = LB(1) - MW
      DO 50 J = 1, LWM
      LB(J+1) = LB(J)/256
      LB(J) = LB(J) - LB(J+1)*256
   50 CONTINUE
      IF (LIN(L).LT.0) LB(LW) = LB(LW) + 128
      DO 60 J = 1, MR
      CHR(I1 + J) = CHAR(LB(LWP - J))
   60 CONTINUE
   70 RETURN
      END
*
*****************************************************
* CONVCI -- CONVERTS CHARACTER DATA TO PACKED INTEGER
*****************************************************
*
      SUBROUTINE CONVCI(CHR,LOUT,LEN,N)
C THIS SUBROUTINE CONVERTS CHARACTER DATA TO PACKED INTEGER
C SHOULD HAVE N = (LEN - 1)/LW + 1
      CHARACTER*1 CHR(*)
      DIMENSION LOUT(N)
      SAVE LW,NW
      DATA LW,NW /4,-2147483647/
      IF (LEN.LT.1) GO TO 60
      MW = NW - 1
      L = (LEN - 1)/LW + 1
      M = LEN/LW
      MR = LEN - M*LW
      LW1 = LW - 1
      IF (M.EQ.0) GO TO 30
      DO 20 J = 1, M
      J1 = (J - 1)*LW + 1
      ITC = ICHAR(CHR(J1))
      LOUT(J) = ITC
      IF (ITC.GE.128) LOUT(J) = LOUT(J) - 128
      DO 10 I = 1, LW1
      LOUT(J) = ICHAR(CHR(J1+I)) + 256*LOUT(J)
   10 CONTINUE
      IF (ITC.GE.128) LOUT(J) = LOUT(J) + MW
   20 CONTINUE
   30 IF (MR.EQ.0) GO TO 60
      J1 = M*LW + 1
      ITC = ICHAR(CHR(J1))
      LOUT(L) = ITC
      IF (ITC.GE.128) LOUT(L) = LOUT(L) - 128
      IF (MR.EQ.1) GO TO 50
      MR1 = MR - 1
      DO 40 I = 1, MR1
      LOUT(L) = ICHAR(CHR(J1+I)) + 256*LOUT(L)
   40 CONTINUE
   50 IT1 = 256**(LW*L - LEN)
      LOUT(L) = IT1*LOUT(L)
      IF (ITC.GE.128) LOUT(L) = LOUT(L) + MW
   60 RETURN
      END
*
*****************************************************
* ENCODE -- ENCODES BINARY TO ASCII
*****************************************************
*
      SUBROUTINE ENCODE (LIN,LEN)
C THIS SUBROUTINE ENCODES BINARY TO ASCII
C WRITTEN FOR THE IBM 3090 VF - VIKTOR K. DECYK, UCLA
C DIMENSION LB(LW), IA(IB+1), LA > = 64
C MW = -2**(8*LW-1)
      DIMENSION LIN(1)
      DIMENSION LB(4), IA(4)
      SAVE LW,NW,IB,LA,IS,IA,K
      DATA LW,NW,IB,LA,IS,K /4,-2147483647,3,83,34,1/
      IF (LEN.LT.1) GO TO 70
      MW = NW - 1
      L = (LEN - 1)/LW + 1
      LWM = LW - 1
      LWP = LW + 1
      IC = IB + 1
      LS = 256**(LW - IC)
      DO 60 I = 1, L
      LB(1) = LIN(I)
      IF (LB(1).LT.0) LB(1) = LB(1) - MW
      DO 10 J = 1, LWM
      LB(J+1) = LB(J)/256
      LB(J) = LB(J) - LB(J+1)*256
   10 CONTINUE
      IF (LIN(I).LT.0) LB(LW) = LB(LW) + 128
      DO 50 J = 1, LW
      IA(K) = LB(LWP - J)
      II = IA(K)
C     IF (II.EQ.0) GO TO 50
      K = K + 1
      IF (K.LE.IB) GO TO 50
      LINE = 0
      DO 20 JJ = 1, IB
      LINE = IA(JJ) + 256*LINE
   20 CONTINUE
      IA(IC) = LINE
      DO 30 JJ = 1, IB
      J1 = IC - JJ
      IA(J1) = IA(J1+1)/LA
      IA(J1+1) = IA(J1+1) - IA(J1)*LA
   30 CONTINUE
      LINE = 0
      DO 40 JJ = 1, IC
      IA(JJ) = IA(JJ) + IS
      IF (IA(JJ).GT.90) IA(JJ) = IA(JJ) + 6
      LINE = IA(JJ) + 256*LINE
   40 CONTINUE
      LINE = LS*LINE
      CALL BUFFPK (LINE,IC)
      K = 1
   50 CONTINUE
   60 CONTINUE
      GO TO 120
   70 IC = K - 1
      IF (IC.LT.1) GO TO 110
      LINE = 0
      DO 80 JJ = 1, IC
      LINE = IA(JJ) + 256*LINE
   80 CONTINUE
      IC1 = IC + 1
      IA(IC1) = LINE
      DO 90 JJ = 1, IC
      J1 = IC1 - JJ
      IA(J1) = IA(J1+1)/LA
      IA(J1+1) = IA(J1+1) - IA(J1)*LA
   90 CONTINUE
      LINE = 0
      DO 100 JJ = 1, IC1
      IA(JJ) = IA(JJ) + IS
      IF (IA(JJ).GT.90) IA(JJ) = IA(JJ) + 6
      LINE = IA(JJ) + 256*LINE
  100 CONTINUE
      LINE = LINE*256**(LW - IC1)
      CALL BUFFPK (LINE,IC1)
      K = 1
      IC = 0
  110 CALL BUFFPK (LINE,IC)
  120 RETURN
      END
*
*****************************************************
* BUFFPK
*****************************************************
*
      SUBROUTINE BUFFPK(LINE,N)
C DIMENSION LOUT((LMAX-1)/LW+2)
      DIMENSION LINE(1)
      DIMENSION LOUT(19)
      SAVE LEN,LMAX,LW,LOUT
      DATA LEN,LMAX,LW /0,72,4/
      DATA LOUT(1) /0/
      NC = N
      NCR = N
      IF (N.GT.0) GO TO 10
      IF (LEN.GT.0) GO TO 50
      GO TO 70
   10 NCR = (NC + LEN) - LMAX
      IF (NCR.GE.0) NC = LMAX - LEN
      M = LEN/LW
      L = LEN - LW*M
      LR = LW - L
      LS = 256**L
      LL = 1
      IF (L.GT.0) LL = 256**LR
      NCT = NC
      I = 1
   20 LT = LINE(I)/LS
      LOUT(M+I) = LOUT(M+I) + LT
      IF (NCT.LT.LR) GO TO 40
      LOUT(M+I+1) = (LINE(I) - LS*LT)*LL
   30 IF (NCT.LE.LW) GO TO 40
      I = I + 1
      NCT = NCT - LW
      GO TO 20
   40 LEN = LEN + NC
      IF (NCR.LT.0) GO TO 70
   50 CALL TPUTC(LOUT,LEN)
      LEN = 0
      IF (NCR.EQ.0) GO TO 60
      NC = NCR
      NCR = NCR - LMAX
      IF (NCR.GE.0) NC = LMAX
      LOUT(1) = LOUT(M+I+1)
      M = -I
      NCT = NCT + NC
      GO TO 30
   60 LOUT(1) = 0
   70 RETURN
      END
*
*****************************************************
* TPUTC  -- TRANSLATES ASCII TO EBCDIC
*****************************************************
*
      SUBROUTINE TPUTC(LOUT,LEN)
C THIS SUBROUTINE TRANSLATES ASCII TO EBCDIC ACCORDING TO THE
C CONVENTIONS AT CORNELL CNSF'S IBM 3090VF AND WRITES RESULT TO METAFILE
C VIKTOR K. DECYK, UCLA
C DIMENSION LB(LW)
C MW = -2**(8*LW-1)
      DIMENSION LOUT(1)
      DIMENSION IATE(128), LB(4)
      SAVE LW,NW,IATE
   91 FORMAT (18A4)
      DATA LW,NW /4,-2147483647/
C EBCDIC CODE FOR ASCII 124 IS NON-STANDARD
      DATA IATE /0,1,2,3,55,45,46,47,22,5,37,11,12,13,14,15,16,17,18,19,
     160,61,50,38,24,25,63,39,34,29,53,31,64,90,127,123,91,108,80,125,77
     2,93,92,78,107,96,75,97,240,241,242,243,244,245,246,247,248,249,122
     3,94,76,126,110,111,124,193,194,195,196,197,198,199,200,201,209,210
     4,211,212,213,214,215,216,217,226,227,228,229,230,231,232,233,173,2
     524,189,95,109,121,129,130,131,132,133,134,135,136,137,145,146,147,
     6148,149,150,151,152,153,162,163,164,165,166,167,168,169,192,79,208
     7,161,7/
      IF (LEN.LT.1) GO TO 40
      MW = NW - 1
      L = (LEN - 1)/LW + 1
      LW1 = LW - 1
      DO 30 I = 1, L
      LB(1) = LOUT(I)
      DO 10 J = 1, LW1
      LB(J+1) = LB(J)/256
      LB(J) = LB(J) - LB(J+1)*256
   10 CONTINUE
      IT1 = IATE(LB(LW)+1)
      IT2 = IT1
      IF (IT1.GE.128) IT1 = IT1 - 128
      DO 20 J = 1, LW1
      IT1 = IATE(LB(LW-J)+1) + 256*IT1
   20 CONTINUE
      LOUT(I) = IT1
      IF (IT2.GE.128) LOUT(I) = LOUT(I) + MW
   30 CONTINUE
      WRITE (19,91) (LOUT(J),J=1,L)
   40 RETURN
      END
*
*****************************************************
* startg -- initializes compressed raster device
*****************************************************
*
      SUBROUTINE STARTG
C THIS SUBROUTINE INITIALIZES COMPRESSED RASTER DEVICE
      COMMON /GRDEV/ LX, LY, NBIT, IFRMT, INTRL, IXOR
      CHARACTER*1 C
      DIMENSION LXS(7), LYS(7), NBITS(4)
      SAVE LXS,LYS,NBITS,ISTART
   91 FORMAT (61H ENTER (1,2,3,4) FOR (COMPRESSED,ENCODED,BINARY,MFE) FO
     1RMAT: )
   92 FORMAT (A1)
   93 FORMAT (75H ENTER (1,2,3,4,5,6,7) FOR (MACPLUS,MACII,CGA,EGA,3179G
     1,ASCII,OTHER) SIZE: )
   94 FORMAT (33H ENTER IMAGE SIZE AS LX,LY PAIR: )
   95 FORMAT (42H ENTER (1,2,3,4) FOR (2,4,16,256) COLORS: )
   96 FORMAT (39H ENTER (1,2,3) FOR (CGA,EGA,VGA) SIZE: )
   97 FORMAT (18H PROGRAM EXECUTING)
      DATA LXS /512,640,320,640,720,79,1024/
      DATA LYS /342,480,200,350,384,21,781/
      DATA NBITS /1,2,4,8/
      DATA ISTART /0/
      IF (ISTART.NE.0) GO TO 90
      INTRL = 0
      IXOR = 1
      NPAL = 1
C INQUIRE FORMAT TYPE
10    C = '4'
C     WRITE (6,91)
C     READ (5,92,END=15) C
   15 IFRMT = ICHAR(C) - ICHAR('0')
      IF (IFRMT.EQ.0) GO TO 80
      IF ((IFRMT.LT.1).OR.(IFRMT.GT.4)) GO TO 10
      IF (IFRMT.EQ.4) GO TO 60
      IF (IFRMT.NE.2) GO TO 20
C     CLOSE(UNIT=10)
C     OPEN(UNIT=10,FILE='MOVIEF',FORM='FORMATTED',STATUS='UNKNOWN')
C NOT MFE FORMAT
C FIRST, INQUIRE SCREEN SIZE
20    continue
C  20 WRITE (6,93)
      READ (5,92,END=25) C
   25 ID = ICHAR(C) - ICHAR('0')
      IF (ID.EQ.0) GO TO 80
      IF ((ID.LT.1).OR.(ID.GT.7)) GO TO 20
      LX = LXS(ID)
      LY = LYS(ID)
      IF (ID.LT.7) GO TO 40
C VARIABLE SCREEN SIZE
30    continue
C  30 WRITE (6,94)
      READ (5,*,END=35) LX, LY
   35 IF ((LX.LT.1).OR.(LX.GT.LXS(7))) GO TO 30
      IF ((LY.LT.1).OR.(LY.GT.LYS(7))) GO TO 30
C NEXT, INQUIRE NUMBER OF COLOR BITS
40    continue
C  40 WRITE (6,95)
      READ (5,92,END=45) C
   45 ID = ICHAR(C) - ICHAR('0')
      IF (ID.EQ.0) GO TO 80
      IF ((ID.LT.1).OR.(ID.GT.4)) GO TO 40
      NBIT = NBITS(ID)
      GO TO 70
C MFE FORMAT
60    continue
C  60 WRITE (6,96)
      READ (5,92,END=65) C
   65 ID = ICHAR(C) - ICHAR('0')
      IF (ID.EQ.0) GO TO 80
      IF ((ID.LT.1).OR.(ID.GT.3)) GO TO 60
      IF (ID.EQ.1) NBIT = 2
      IF (ID.EQ.2) NBIT = 1
      IF (ID.EQ.3) NBIT = 8
      IF (ID.EQ.1) INTRL = 1
      IF (ID.EQ.3) NPAL = 0
      IF (ID.LT.3) ID = ID + 2
      LX = LXS(ID)
      LY = LYS(ID)
   70 IF ((IFRMT.EQ.1).OR.(IFRMT.EQ.2)) IXOR = 1
      IF ((IFRMT.EQ.3).OR.(IFRMT.EQ.4)) IXOR = 0
      CALL HEADER(IFRMT,LX,LY,NBIT)
      IF (NPAL.EQ.0) CALL WRPAL(C,NPAL,IFRMT)
      ISTART = 1
      GO TO 90
   80 STOP 1
   90 WRITE (6,97)
      END
*
*****************************************************
* GRPARM -- INITIALIZES GRAPHICS SIZE PARAMETERS
*****************************************************
*
      SUBROUTINE GRPARM(IPARM,IRX,IRY,MNX,MNY,LNX,LNY,ICH,ICW,ISLB,ISTCX
     1,ISTCY,NTCX,NTCY,IGSTYL)
C THIS SUBROUTINE INITIALIZES GRAPHICS SIZE PARAMETERS
C IPARM = (1,2) = PARAMETERS FOR (LINE,CONTOUR) PLOTS
      COMMON /GRDEV/ LX, LY, NBIT, IFRMT, INTRL, IXOR
      IF (IPARM.LE.0) GO TO 40
      IRX = LX
      IRY = LY
      XS = FLOAT(LX)/1024.
      YS = FLOAT(LY)/780.
      ICH = 26.*YS + .5
      ICW = 18.*XS + .5
      IGSTYL = 1
      IF (IPARM.GT.1) GO TO 20
C PARAMETERS FOR LINE PLOTS
      MNX = 225.*XS + .5
      MNY = 115.*YS + .5
      LNX = 785.*XS + .5
      LNY = 650.*YS + .5
      ISLB = 225.*XS + .5
      ISTCX = 10.*XS + .5
      ISTCY = 10.*YS + .5
      NTCX = 10
      NTCY = 10
      GO TO 40
   20 IF (IPARM.GT.2) GO TO 30
C PARAMETERS FOR CONTOUR PLOTS
      MNX = 225.*XS + .5
      MNY = 80.*YS + .5
      LNX = 690.*XS + .5
      LNY = 690.*YS + .5
      ISLB = 225.*XS + .5
      ISTCX = 0
      ISTCY = 0
      NTCX = 1
      NTCY = 1
      GO TO 40
   30 IF (IPARM.GT.3) GO TO 40
C PARAMETERS FOR TEXT PLOTS
      MNX = 0
      MNY = 0
      LNX = LX - 1
      LNY = LY - 1
      ISLB = 0
      ISTCX = 0
      ISTCY = 0
      NTCX = 1
      NTCY = 1
*     **************************************************
*     WHAT ABOUT PARAMETERS FOR RASTER IMAGES???
*     WHAT SHOULD THESE BE??????
*     **************************************************
   40 RETURN
      END
*
*****************************************************
* DRAWG -- GENERIC DRAW ROUTINE
*****************************************************
*
*
*****************************************************
* drawg -- generic draw routine
*****************************************************
*
      SUBROUTINE DRAWG(C,X,Y,IC,LWTYPE,L,ICODE)
C THIS SUBROUTINES IS A GENERIC DRAW ROUTINE
C IF ICODE = 0 ZEROES IMAGE
C IF ICODE = 1 PERFORMS MOVE
C IF ICODE = 2 DRAWS LINE
C IF ICODE = 3 DRAWS POINT
C IF ICODE = 4 DRAWS DASHED LINE
C IF ICODE = 5 WRITES CHARACTERS
C IF ICODE = 6 DRAWS IMAGE
C LX = NUMBER OF PIXELS IN X, LY = NUMBER OF PIXELS IN Y
C NBIT = NUMBER OF BITS PER PIXEL, NBIT < 9
      PARAMETER(LXM=1024,LYM=781,NBITD=8)
      PARAMETER(NPIXD=8/NBITD,LZM=(LXM-1)/NPIXD+1,LENBM=LZM*LYM)
      PARAMETER(LENGM=LENBM+LZM,LZGM=LZM+1)
      COMMON /GRDEV/ LX, LY, NBIT, IFRMT, INTRL, IXOR
      COMMON /RCOM/ MINX,LENX,MINY,LENY,ICX,ICY,
     1XMIN,XMAX,YMIN,YMAX,DX,DY,CSX,CSY
      CHARACTER*1 C
      CHARACTER*1 BLANK
      CHARACTER*1 G(LXM*LYM)
      CHARACTER*1 IMAGE(LENBM), JMAGE(LENBM)
      CHARACTER*1 IMG(LENGM), LINE(LZM), LIMG(LZGM)
      CHARACTER*1 CTYPE
      DIMENSION ICOLOR(16), ITWOBC(16)
      DIMENSION IETA(256)
      SAVE ISTART,ICC,G,CTYPE,BLANK,ICOLOR,ITWOBC,IETA
      SAVE LZ,LENB,LENG,LZG,JMAGE
   91 FORMAT (24H SIZE ERROR, LXM, LYM = ,2I12,10H LX, LY = ,2I12)
      DATA ISTART,ICC /0,-1/
C COLOR TABLE FOR 4 BIT COLOR
      DATA ICOLOR /0,1,2,3,4,5,6,7,1,1,2,3,4,5,6,7/
C COLOR TABLE FOR 2 BIT COLOR
      DATA ITWOBC /0,1,1,1,2,2,3,3,1,1,1,1,2,2,3,3/
C EBCDIC/ASCII TRANSLATION WITH CONVENTIONS AT CORNELL CNSF IBM 3090VF.
C ASCII CODES FOR EBCDIC 74,79,113,139,155 ARE NON-STANDARD
C EBCDIC CODES 28,30,106 ARE ADDED FOR IBM COMPATIBILITY
      DATA IETA /0,1,2,3,-1,9,-1,127,-1,-1,-1,11,12,13,14,15,16,17,18,19
     1,-1,-1,8,-1,24,25,-1,-1,28,29,30,31,-1,-1,28,-1,-1,10,23,27,-1,-1,
     2-1,-1,-1,5,6,7,-1,-1,22,-1,-1,30,-1,4,-1,-1,-1,-1,20,21,-1,26,32,-
     31,-1,-1,-1,-1,-1,-1,-1,-1,92,46,60,40,43,124,38,-1,-1,-1,-1,-1,-1,
     4-1,-1,-1,33,36,42,41,59,94,45,47,-1,-1,-1,-1,-1,-1,-1,-1,124,44,37
     5,95,62,63,-1,94,-1,-1,-1,-1,-1,-1,-1,96,58,35,64,39,61,34,-1,97,98
     6,99,100,101,102,103,104,105,-1,123,-1,-1,-1,-1,-1,106,107,108,109,
     7110,111,112,113,114,-1,125,-1,-1,-1,-1,-1,126,115,116,117,118,119,
     8120,121,122,-1,-1,-1,91,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
     9-1,93,-1,-1,123,65,66,67,68,69,70,71,72,73,-1,-1,-1,-1,-1,-1,125,7
     A4,75,76,77,78,79,80,81,82,-1,-1,-1,-1,-1,-1,92,-1,83,84,85,86,87,8
     B8,89,90,-1,-1,-1,-1,-1,-1,48,49,50,51,52,53,54,55,56,57,-1,-1,-1,-
     C1,-1,-1/
      DATA CTYPE /' '/
      IF (ISTART.NE.0) GO TO 30
C INITIALIZE VARIABLES
      IF ((LX.LE.LXM).AND.(LY.LE.LYM)) GO TO 10
      WRITE (6,91) LXM, LYM, LX, LY
      STOP 1
   10 NPIX = 8/NBIT
      LZ = (LX - 1)/NPIX + 1
      LENB = LZ*LY
      LENG = LENB + LZ
      LZG = LZ + 1
      BLANK = CHAR(0)
      DO 20 J = 1, LENB
      JMAGE(J) = BLANK
   20 CONTINUE
      ISTART = 1
C CHECK CODE
   30 IF ((ICODE.LT.0).OR.(ICODE.GT.6)) GO TO  80
      IF (ICODE.GT.0) GO TO 40
C CLEAR IMAGE
      CALL ZIMAGE(G,BLANK,LX,LY)
      GO TO  80
   40 IF (ICODE.GT.5) GO TO 70
C PERFORM DRAW INSTRUCTION
      IF (IC.EQ.ICC) GO TO 50
      ICC = IC
      ICT = IC - (IC/16)*16
      IF ((NBIT.EQ.1).AND.(ICT.GT.0)) ICT = 1
      IF (NBIT.EQ.2) ICT = ITWOBC(ICT+1)
      IF (NBIT.GE.4) ICT = ICOLOR(ICT+1)
      CTYPE = CHAR(ICT)
C FIRST PERFORM SCALING TO RASTER UNITS AND CLIP IF NECESSARY
   50 I = (X - XMIN)*DX + .5
      J = (Y - YMIN)*DY + .5
      IF (I.LT.0) I = 0
      IF (I.GT.LENX) I = LENX
      IF (J.LT.0) J = 0
      IF (J.GT.LENY) J = LENY
      I = I + MINX
      J = J + MINY
C PERFORM DRAW
      IF (ICODE.EQ.5) GO TO 60
      IF (ICODE.EQ.1) CALL DMOVE(I,J,ICX,ICY)
      IF (ICODE.EQ.2) CALL DLINE(G,CTYPE,BLANK,LX,LY,I,J,ICX,ICY,LWTYPE)
      IF (ICODE.EQ.3) CALL DPNT(G,CTYPE,BLANK,LX,LY,I,J,ICX,ICY,LWTYPE)
      IF (ICODE.EQ.4) CALL DASHLN(G,CTYPE,BLANK,LX,LY,I,J,ICX,ICY,L,LWTY
     1PE)
      GO TO  80
   60 ICT = IETA(ICHAR(C)+1)
      CALL CDRAW(G,CTYPE,BLANK,LX,LY,ICT,I,J,ICX,ICY,CSX,CSY,LWTYPE)
      GO TO  80
   70 CALL GIMAGE (G,IMAGE,LX,LY,LZ,NBIT,INTRL)
      CALL PCTSAV(IMAGE,JMAGE,IMG,LINE,LIMG,IXOR,LZ,LY,LENB,LENG,LZG,IFR
     1MT,INTRL)
   80 RETURN
      END
*
*****************************************************
* COMPRG -- COMPRESSES RASTER IMAGE AND WRITES TO DISK
*****************************************************
*
      SUBROUTINE COMPRG( G, NX, NY, NXV, inputbits, filter, irev )

C THIS SUBROUTINE COMPRESSES RASTER IMAGE AND WRITES RESULT TO DISK
C INPUT IS IN ARRAY G, AND OUTPUT IS WRITTEN TO DISK
C IF LX < NX OR LY < NY, THEN IMAGE IS TRUNCATED
C IF LX > NX OR LY > NY, THEN IMAGE IS PADDED WITH NULLS
C LX = NUMBER OF PIXELS IN X, LY = NUMBER OF PIXELS IN Y
C NBIT = NUMBER OF BITS PER PIXEL, NBIT < 9

      PARAMETER(LXM=1024,LYM=781,NBITD=8)
      PARAMETER(NPIXD=8/NBITD,LZM=(LXM-1)/NPIXD+1,LENBM=LZM*LYM)
      PARAMETER(LENGM=LENBM+LZM,LZGM=LZM+1)
      COMMON /GRDEV/ LX, LY, NBIT, IFRMT, INTRL, IXOR
      CHARACTER*1 G(NXV,NY)
      integer     inputbits
      integer     filter
      CHARACTER*1 BLANK
      CHARACTER*1 IMAGE(LENBM), JMAGE(LENBM)
      CHARACTER*1 IMG(LENGM), LINE(LZM), LIMG(LZGM)
      SAVE ISTART,BLANK
      SAVE LZ,LENB,LENG,LZG,JMAGE
   91 FORMAT (24H SIZE ERROR, LXM, LYM = ,2I12,10H LX, LY = ,2I12)
      DATA ISTART /0/

      IF (ISTART.NE.0) GO TO 30

C INITIALIZE VARIABLES
      IF ((LX.LE.LXM).AND.(LY.LE.LYM)) GO TO 10
      WRITE (6,91) LXM, LYM, LX, LY
      STOP 1
   10 NPIX = 8/NBIT
      LZ = (LX - 1)/NPIX + 1
      LENB = LZ*LY
      LENG = LENB + LZ
      LZG = LZ + 1
      BLANK = CHAR(0)
         DO 20 J = 1, LENB
            JMAGE(J) = BLANK
   20    CONTINUE
      ISTART = 1

*     *****************************************************
*     apply filter if needed
*        there are 4 filters available:
*           filter = 0  -- do not filter
*           filter = 1  -- all non-background colors to white
*           filter = 2  -- threshold
*           filter = 3  -- dither
*        there are 3 output formats available:
*           nbit = 1    -- monocrome
*           nbit = 2    -- 4 color
*           nbit = 8    -- 256 color
*     *****************************************************

   30 continue

*     *****************************************************
*     if filtering not requested do nothing
*     *****************************************************

      if( filter .eq. 0 ) then

*     *****************************************************
*     if the number of bits in the output image is greate
*     than or equal to those in the input image, do nothing
*     *****************************************************

      else if( nbit .ge. inputbits ) then

*     *****************************************************
*     filtering is all non-background colors to white
*     *****************************************************

      else if( filter .eq. 1 ) then
         call towhite( g, nx, ny, inputbits )

*     *****************************************************
*     filtering is by threshold
*     *****************************************************

      else if( filter .eq. 2 ) then
         if( nbit .eq. 1 ) then
            call thrhld1( g, nx, ny, inputbits )
         else if( nbit .eq. 2 ) then
            call thrhld2( g, nx, ny )
         else
         end if

*     *****************************************************
*     filtering is by dither
*     *****************************************************

      else
         if( nbit .eq. 1 ) then
            if( inputbits .eq. 2 ) then
               call dithr12( g, nx, ny )
            else
               call dithr18( g, nx, ny )
            end if
         else if( nbit .eq. 2 ) then
            call dithr28( g, nx, ny )
         else
         end if
      end if

*     *****************************************************
*     now compress and write out the image
*     *****************************************************

      CALL GRIMAGE (G,IMAGE,NX,NY,NXV,LX,LY,LZ,NBIT,INTRL,irev)
      CALL PCTSAV(IMAGE,JMAGE,IMG,LINE,LIMG,IXOR,LZ,LY,LENB,LENG,LZG,IFR
     1MT,INTRL)
      RETURN
      END
*
*****************************************************
* towhite -- turns all non-background colors to white
*****************************************************
*
      subroutine towhite( g, nx, ny, inputbits )

      COMMON /GRDEV/ LX, LY, NBIT, IFRMT, INTRL, IXOR

      character*1 g(nx,ny)


*        ********************************
*        convert each byte in the g array
*        ********************************

         do 1 ix = 1,nx
            do 2 iy = 1,ny

*              *****************************************
*              since we are only reducing the # of bits,
*              then if the inputbits = 2, the output
*              bits must = 1.
*              if inputbits = 2, there are 4 input colors
*              with 0 = background
*              for bits = 1, 0 is background, 1 = the color
*              ******************************************

               if( inputbits .eq. 2 ) then
                  is = ichar( g( ix,iy ) )
                  g( ix,iy ) = char( is - (is/4)*4 )
                  if( g( ix,iy ) .ne. char( 0 ) )
     .               g( ix,iy ) = char( 1 )

*              **************************************************
*              if inputbits = 8, there are 256 input colors
*              with 0 = background, all others are color
*              if nbit is 1, the non-background color is 1
*              if nbit is 2, the non-background color is 3
*              **************************************************

               else
                  if( nbit .eq. 1 ) then
                     if( g( ix,iy ) .ne. char( 0 ) )
     .                  g( ix,iy ) = char( 1 )
                  else
                     if( g( ix,iy ) .ne. char( 0 ) )
     .                  g( ix,iy ) = char( 3 )
                  end if

               end if

2           continue
1        continue
      return
      end
*
*****************************************************
* thrhld1 -- threshold to 1 bit output
*****************************************************
*
      subroutine thrhld1( g, nx, ny, inputbits )

      character*1 g(nx,ny)


*     **************************************************
*     find the max and min values in g for 8 bit images
*     **************************************************

      if( inputbits .eq. 8 ) then
         min = 255
         max = 0

            do 3 ix = 1,nx
               do 4 iy = 1,ny
                  ival = ichar( g(ix,iy) )
                  if( ival .lt. min )
     .               min = ival
                  if( ival .gt. max )
     .               max = ival
4              continue
3           continue

         idist = max - min + 1
         ihalf = idist / 2
         midpt = min + ihalf - 1

      end if

*        ********************************
*        convert each byte in the g array
*        ********************************

         do 1 ix = 1,nx
            do 2 iy = 1,ny

*              *****************************************
*              if inputbits = 2, convert 0 and 1 to 0
*              and 2 and 3 to 1
*              ******************************************

               if( inputbits .eq. 2 ) then

*                 **************************************
*                 only except stuff in the low order 2
*                 bits of g
*                 ***************************************

                  is = ichar( g( ix,iy ) )
                  g( ix,iy ) = char( is - (is/4)*4 )

*                 **************************************
*                 now compare with the threshold
*                 ***************************************

                  if( g( ix,iy ) .le. char( 1 ) ) then
                     g( ix,iy ) = char( 0 )
                  else
                     g( ix,iy ) = char( 1 )
                  end if

*              **************************************************
*              if inputbits = 8, there are 256 input colors
*              convert 0 - midpt to 0 and midpt - 255 to 1
*              **************************************************

               else
                  if( ichar( g(ix,iy) ) .le. midpt ) then
                     g( ix,iy ) = char( 0 )
                  else
                     g( ix,iy ) = char( 1 )
                  end if

               end if

2           continue
1        continue
      return
      end
*
*****************************************************
* thrhld2 -- threshold to 2 bit output
*            2 bit output has the following colors:
*               0 = background
*               1 = magenta
*               2 = cyan
*               3 = white
*****************************************************
*
      subroutine thrhld2( g, nx, ny )

      character*1 g(nx,ny)
      integer     midpt( 3 )


*     ***********************************************
*     we can only be dealing with 8 bit input that
*     we are taking down to 2 bit output
*     ***********************************************

*     **************************************************
*     find the max and min values in g
*     **************************************************

      min = 255
      max = 0

         do 3 ix = 1,nx
            do 4 iy = 1,ny
               ival = ichar( g(ix,iy) )
               if( ival .lt. min )
     .            min = ival
               if( ival .gt. max )
     .            max = ival
4           continue
3        continue

      idist = max - min + 1
      ihalf = idist / 2
      midpt( 2 ) = min + ihalf - 1
      ihalf = idist / 4
      midpt( 1 ) = min + ihalf - 1
      midpt( 3 ) = midpt( 2 ) + ihalf - 1

*        ********************************
*        convert each byte in the g array
*        ********************************

         do 1 ix = 1,nx
            do 2 iy = 1,ny

*              **************************************************
*              if inputbits = 8, there are 256 input colors
*              convert 0 - 63 to 0, 64 - 127 to 1
*                      128 - 191 to 2, 192 - 255 to 3
*              **************************************************

               if( ichar( g(ix,iy) ) .le. midpt( 1 ) ) then
                  g( ix,iy ) = char( 0 )
               else if( ichar( g(ix,iy) ) .le. midpt( 2 ) ) then
                  g( ix,iy ) = char( 1 )
               else if( ichar( g(ix,iy) ) .le. midpt( 3 ) ) then
                  g( ix,iy ) = char( 2 )
               else
                  g( ix,iy ) = char( 3 )
               end if

2           continue
1        continue
      return
      end
*
*****************************************************
* dithr12 -- dither a 2 bit input to a 1 bit output
*****************************************************
*
      subroutine dithr12( g, nx, ny )

      character*1 g( nx,ny )
      integer     d( 2,2 )
      data d  / 0, 2,
     .          3, 1  /


*        ********************************
*        convert each byte in the g array
*        ********************************

         do 1 ix = 1,nx
            ixd = mod( ix, 2 ) + 1

               do 2 iy = 1,ny
                  iyd = mod( iy, 2 ) + 1

*                 **************************************
*                 only except stuff in the low order 2
*                 bits of g
*                 ***************************************

                  is = ichar( g( ix,iy ) )
                  g( ix,iy ) = char( is - (is/4)*4 )

*                 **************************************
*                 now compare with the dither threshold
*                 ***************************************

                  if( g( ix, iy ) .le. char( d( ixd,iyd ) ) ) then
                     g( ix, iy ) = char( 0 )
                  else
                     g( ix, iy ) = char( 1 )
                  end if

2              continue
1        continue
      return
      end
*
*****************************************************
* dithr18 -- dither an 8 bit input to a 1 bit output
*****************************************************
*
      subroutine dithr18( g, nx, ny )

      character*1 g( nx,ny )
      integer     d( 8,8 )
      data d  /  0, 32,  8, 40,  2, 34, 10, 42,
     .          48, 16, 56, 24, 50, 18, 58, 26,
     .          12, 44,  4, 36, 14, 46,  6, 38,
     .          60, 28, 52, 20, 62, 30, 54, 22,
     .           3, 35, 11, 43,  1, 33,  9, 41,
     .          51, 19, 59, 27, 49, 17, 57, 25,
     .          15, 47,  7, 39, 13, 45,  5, 37,
     .          63, 31, 55, 23, 61, 29, 53, 21   /


*     **************************************************
*     find the max and min values in g
*     **************************************************

      min = 255
      max = 0

         do 3 ix = 1,nx
            do 4 iy = 1,ny
               ival = ichar( g(ix,iy) )
               if( ival .lt. min )
     .            min = ival
               if( ival .gt. max )
     .            max = ival
4           continue
3        continue

      if( min .lt. 64 ) then
         min = 0
      else if( min .lt. 128 ) then
         min = 64
      else if( min .lt. 192 ) then
         min = 128
      else
         min = 192
      end if

      if( max .lt. 63 ) then
         max = 63
      else if( max .lt. 127 ) then
         max = 127
      else if( max .lt. 191 ) then
         max = 191
      else
         max = 255
      end if

      if( max - min .eq. 63 ) then
         mfactor = 1
      else if( max - min .eq. 127 ) then
         mfactor = 2
      else if( max - min .eq. 191 ) then
         mfactor = 3
      else
         mfactor = 4
      end if

*        ********************************
*        convert each byte in the g array
*        ********************************

         do 1 ix = 1,nx
            ixd = mod( ix, 8 ) + 1

               do 2 iy = 1,ny
                  iyd = mod( iy, 8 ) + 1

                  if(  ichar( g(ix,iy) )
     .                 .le.
     .                 ( d(ixd,iyd) * mfactor ) + min  ) then
                     g( ix, iy ) = char( 0 )
                  else
                     g( ix, iy ) = char( 1 )
                  end if

2              continue
1        continue
      return
      end
*
*****************************************************
* dithr28 -- dither an 8 bit input to a 2 bit output
*****************************************************
*
      subroutine dithr28( g, nx, ny )

      common /dithpal/ pal64, npal64
      character*1 pal64( 3,256 )

      character*1 g( nx,ny )

      character*1 red( 1000,1000 )
      character*1 green( 1000,1000 )
      character*1 blue( 1000,1000 )

      integer     d( 8,8 )
      data d  /  0, 32,  8, 40,  2, 34, 10, 42,
     .          48, 16, 56, 24, 50, 18, 58, 26,
     .          12, 44,  4, 36, 14, 46,  6, 38,
     .          60, 28, 52, 20, 62, 30, 54, 22,
     .           3, 35, 11, 43,  1, 33,  9, 41,
     .          51, 19, 59, 27, 49, 17, 57, 25,
     .          15, 47,  7, 39, 13, 45,  5, 37,
     .          63, 31, 55, 23, 61, 29, 53, 21   /


*     ********************************************************
*     if there is no palette present, use the default palette
*     ********************************************************


      if( npal64 .eq. 0 ) then
         pal64( 1, 1 ) = char(  0 )
         pal64( 2, 1 ) = char(  0 )
         pal64( 3, 1 ) = char(  0 )

         pal64( 1, 2 ) = char(  0 )
         pal64( 2, 2 ) = char(  0 )
         pal64( 3, 2 ) = char( 63 )

         pal64( 1, 3 ) = char(  0 )
         pal64( 2, 3 ) = char( 63 )
         pal64( 3, 3 ) = char(  0 )

         pal64( 1, 4 ) = char(  0 )
         pal64( 2, 4 ) = char( 63 )
         pal64( 3, 4 ) = char( 63 )

         pal64( 1, 5 ) = char( 63 )
         pal64( 2, 5 ) = char(  0 )
         pal64( 3, 5 ) = char(  0 )

         pal64( 1, 6 ) = char( 63 )
         pal64( 2, 6 ) = char(  0 )
         pal64( 3, 6 ) = char( 63 )

         pal64( 1, 7 ) = char( 63 )
         pal64( 2, 7 ) = char( 63 )
         pal64( 3, 7 ) = char(  0 )

         pal64( 1, 8 ) = char( 63 )
         pal64( 2, 8 ) = char( 63 )
         pal64( 3, 8 ) = char( 63 )

         npal64 = 8

      end if

*     **************************************
*     do this if the palet has > 4 entries
*     **************************************


      if( npal64 .gt. 4 ) then

*        ***********************************************
*        look up each pixel in the image on the palette
*        and create a separate red, green, and blue image
*        ***********************************************

         do 1 ix = 1,nx
            do 2 iy = 1,ny

               ipalno = ichar( g( ix,iy ) ) + 1
               if( ipalno .gt. npal64 ) then
                  red(   ix,iy ) = char( 0 )
                  green( ix,iy ) = char( 0 )
                  blue(  ix,iy ) = char( 0 )
               else
                  red(   ix,iy ) = pal64( 1, ipalno )
                  green( ix,iy ) = pal64( 2, ipalno )
                  blue(  ix,iy ) = pal64( 3, ipalno )
               end if
2           continue
1        continue

*        ***********************************************
*        dither the reds
*        ***********************************************

         do 11 ix = 1,nx
            ixd = mod( ix, 8 ) + 1

               do 12 iy = 1,ny
                  iyd = mod( iy, 8 ) + 1

                  if( red( ix, iy ) .le. char( d( ixd, iyd ) ) ) then
                     red( ix, iy ) = char( 0 )
                  else
                     red( ix, iy ) = char( 1 )
                  end if

12             continue
11       continue

*        ***********************************************
*        dither the greens
*        ***********************************************

         do 21 ix = 1,nx
            ixd = mod( ix, 8 ) + 1

               do 22 iy = 1,ny
                  iyd = mod( iy, 8 ) + 1

                  if( green( ix, iy ) .le. char( d( ixd, iyd ) ) ) then
                     green( ix, iy ) = char( 0 )
                  else
                     green( ix, iy ) = char( 1 )
                  end if

22             continue
21       continue

*        ***********************************************
*        dither the blues
*        ***********************************************

         do 31 ix = 1,nx
            ixd = mod( ix, 8 ) + 1

               do 32 iy = 1,ny
                  iyd = mod( iy, 8 ) + 1

                  if( blue( ix, iy ) .le. char( d( ixd, iyd ) ) ) then
                     blue( ix, iy ) = char( 0 )
                  else
                     blue( ix, iy ) = char( 1 )
                  end if

32             continue
31       continue

*        ***********************************************
*        combine the red, green, and blue images back
*        into a single image using the following rules:
*
*           red                    -> magenta = 1
*           green                  -> cyan = 2
*           blue                   -> cyan = 2
*           red + green            -> white = 3
*           red + blue             -> magenta = 1
*           green + blue           -> cyan = 2
*           red + green + blue     -> white = 3
*           nothing                -> background = 0
*        ***********************************************

         do 5 ix = 1,nx
            do 6 iy = 1,ny

               if( red( ix,iy ) .eq. char( 1 ) ) then
                  if( green( ix,iy ) .eq. char( 1 ) ) then
                     g( ix,iy ) = char( 3 )
                  else
                     g( ix,iy ) = char( 1 )
                  end if
               else if( green( ix,iy ) .eq. char( 1 ) ) then
                  g( ix,iy ) = char( 2 )
               else if( blue( ix,iy ) .eq. char( 1 ) )then
                  g( ix,iy ) = char( 2 )
               else
                  g( ix,iy ) = char( 0 )
               end if

6           continue
5        continue

*     **************************************************
*     if the palette has less then or equal to 4 entries
*     just set any image pixel greater than 3 to 3
*     **************************************************

      else

         do 8 ix = 1,nx
            do 9 iy = 1,ny
               if( g( ix,iy ) .gt. char( 3 ) )
     .            g( ix,iy ) = char( 3 )
9           continue
8        continue

      end if

      return
      end
*
*****************************************************
* CINPUT
*****************************************************
*
      SUBROUTINE CINPUT(C)
C DUMMY CHARACTER READ
      CHARACTER*8 C
      C = '        '
      RETURN
      END
*
*****************************************************
* QUITG -- TERMINATES COMPRESSED RASTER DEVICE
*****************************************************
*
      SUBROUTINE QUITG
C THIS SUBROUTINE TERMINATES COMPRESSED RASTER DEVICE
      COMMON /GRDEV/ LX, LY, NBIT, IFRMT, INTRL, IXOR
      CHARACTER*1 BLANK
      CHARACTER*1 CHR(4)
   91 FORMAT (13H PROGRAM DONE)
      BLANK = CHAR(0)
         DO 10 I = 1, 4
            CHR(I) = BLANK
   10    CONTINUE
      CALL BUFFWR(CHR,4,IFRMT)
      CALL BUFFWR(CHR,0,IFRMT)
      WRITE (6,91)
      RETURN
      END
*
*****************************************************
* VPARSE
*****************************************************
*
C PARSING LIBRARY FOR BEPS1
C COPYRIGHT 1990, REGENTS OF THE UNIVERSITY OF CALIFORNIA
C UPDATE: OCTOBER 16, 1990
      SUBROUTINE VPARSE (CODE,ICODE,CP,IP,AP,INPUT,NC,NCC,NCI,NCR,IVAR,I
     1RC)
C THIS SUBROUTINE PARSES INPUT STRING AND STORES APPROPRIATE VALUE INTO
C APPROPRIATE VARIABLE.  THE VARIABLE NAMES (NC OF THEM) IN THE SYMBOL
C TABLE CODE MUST BE SIX CHARACTERS OR LESS IN LENGTH, AND CREATED IN
C ORDER OF TYPE CHARACTER*8 (NCC OF THEM), INTEGER (NCI OF THEM), AND
C REAL*4 (NCR OF THEM).  INTEGER CODES CORRESPONDING TO THE CHARACTER
C TABLE MUST HAVE BEEN CALCULATED BY CALLS TO SUBROUTINE FINDCC:
C     NC = NCC + NCI + NCR
C     DO 10 I = 1, NC
C     CALL FINDCC(CODE(I),ICODE(I),IC)
C  10 CONTINUE
C FIRST THE ENTRY IVAR IN THE TABLE ICODE IS IDENTIFIED (IRC=1 MEANS
C SYMBOL NOT FOUND).  THEN THE NUMERICAL VALUE OF THE CHARACTERS TO THE
C RIGHT OF THE '=' SIGN IS FOUND (IRC=2 MEANS NO VALID VALUE FOUND), AND
C STORED IN CP, IP, OR AP AS APPROPRIATE FOR CHARACTER, INTEGER, OR REAL
C VARIABLES, RESPECTIVELY.
C     CHARACTER*8 CP(NCC)
C     DIMENSION IP(NCI), AP(NCR)
      CHARACTER*(*) INPUT
      CHARACTER*6 CODE(NC)
      CHARACTER*8 CP(1)
      DIMENSION ICODE(NC)
      DIMENSION IP(1), AP(1)
   91 FORMAT (1X,A6,3H = ,A8)
   92 FORMAT (1X,A6,3H = ,I8)
   93 FORMAT (1X,A6,3H = ,F14.7)
      IM = LEN(INPUT)
      IRC = 1
C FIND THE EQUAL SIGN
      I = 1
   10 IF (INPUT(I:I).EQ.'=') GO TO 20
      I = I + 1
      IF (I.LE.IM) GO TO 10
   20 I = I - 1
C NO VARIABLE LEFT OF EQUAL SIGN
      IF (I.LT.1) GO TO 80
C FIND NUMERICAL CODE FOR VARIABLE NAME
      CALL FINDCC(INPUT(1:I),NUM,IC)
C INVALID CHARACTERS IN NAME
      IF (IC.EQ.0) GO TO 80
C FIND VARIABLE NAME IN TABLE
      J = 0
   30 J = J + 1
      IF (ICODE(J).EQ.NUM) GO TO 40
      IF (J.LT.NC) GO TO 30
C VARIABLE NAME NOT FOUND IN TABLE
      GO TO 80
   40 IVAR = J
      IRC = 2
      IF (I.LT.(IM - 1)) GO TO 50
C NO VALUE RIGHT OF EQUAL SIGN
      GO TO 80
C FIND VARIABLE TYPE
   50 IF (J.GT.NCC) GO TO 60
C CHARACTER VARIABLE
      CP(J) = INPUT(I+2:IM)
      WRITE (6,91) CODE(J), CP(J)
      IRC = 0
      GO TO 80
C FIND VALUE
   60 CALL EVALC(INPUT(I+2:IM),IVAL,VAL,ID)
C NOT A VALID NUMBER
      IF (ID.EQ.0) GO TO 80
C VALID NUMBER
      IRC = 0
      NCT = NCC + NCI
      IF (J.GT.NCT) GO TO 70
C INTEGER VARIABLE
      IP(J-NCC) = IVAL
      WRITE (6,92) CODE(J), IVAL
      GO TO 80
C REAL VARIABLE
   70 AP(J-NCT) = VAL
      WRITE (6,93) CODE(J), VAL
   80 RETURN
      END
*
*****************************************************
* FINDCC
*****************************************************
*
      SUBROUTINE FINDCC(CVAR,NUM,IC)
C THIS SUBROUTINE FINDS NUMERICAL CODE FOR VARIABLE NAME, WHICH CAN
C CONSIST OF LETTERS AND NUMBERS.  UPPER AND LOWER CASE ARE TREATED AS
C EQUIVALENT, AND SPACES ARE IGNORED.  IC IS NUMBER OF SYMBOLS FOUND IN
C VARIABLE NAME (IC = 0 IS RETURNED IF ILLEGAL CHARACTER IS FOUND).
C VARIABLE NAME MUST BE SIX CHARACTERS OR LESS IN LENGTH.
      CHARACTER*1 V
      CHARACTER*(*) CVAR
      DIMENSION IETA(256)
      DATA IB /36/
C ASCII CODES FOR EBCDIC 74,79,113,139,155 ARE NON-STANDARD
C EBCDIC CODES 28,30,106 ARE ADDED FOR IBM COMPATIBILITY
      DATA IETA /0,1,2,3,-1,9,-1,127,-1,-1,-1,11,12,13,14,15,16,17,18,19
     1,-1,-1,8,-1,24,25,-1,-1,28,29,30,31,-1,-1,28,-1,-1,10,23,27,-1,-1,
     2-1,-1,-1,5,6,7,-1,-1,22,-1,-1,30,-1,4,-1,-1,-1,-1,20,21,-1,26,32,-
     31,-1,-1,-1,-1,-1,-1,-1,-1,92,46,60,40,43,124,38,-1,-1,-1,-1,-1,-1,
     4-1,-1,-1,33,36,42,41,59,94,45,47,-1,-1,-1,-1,-1,-1,-1,-1,124,44,37
     5,95,62,63,-1,94,-1,-1,-1,-1,-1,-1,-1,96,58,35,64,39,61,34,-1,97,98
     6,99,100,101,102,103,104,105,-1,123,-1,-1,-1,-1,-1,106,107,108,109,
     7110,111,112,113,114,-1,125,-1,-1,-1,-1,-1,126,115,116,117,118,119,
     8120,121,122,-1,-1,-1,91,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
     9-1,93,-1,-1,123,65,66,67,68,69,70,71,72,73,-1,-1,-1,-1,-1,-1,125,7
     A4,75,76,77,78,79,80,81,82,-1,-1,-1,-1,-1,-1,92,-1,83,84,85,86,87,8
     B8,89,90,-1,-1,-1,-1,-1,-1,48,49,50,51,52,53,54,55,56,57,-1,-1,-1,-
     C1,-1,-1/
      IM = LEN(CVAR)
      NUM = 0
      IC = 0
      I = 0
   10 I = I + 1
      IF ((I.GT.IM).OR.(IC.GE.6)) GO TO 60
      V = CVAR(I:I)
      IF (V.EQ.' ') GO TO 10
C     IV = ICHAR(V) - 64
      IV = IETA(ICHAR(V)+1) - 64
      IF (IV.GT.26) GO TO 30
      IF (IV.LT.1) GO TO 40
   20 NUM = IV + IB*NUM
      IC = IC + 1
      GO TO 10
   30 IV = IV - 32
      IF ((IV.LT.1).OR.(IV.GT.26)) GO TO 50
      GO TO 20
   40 IV = IV + 17
      IF ((IV.LT.1).OR.(IV.GT.10)) GO TO 50
      IV = IV + 26
      GO TO 20
   50 IC = 0
   60 RETURN
      END
*
*****************************************************
* EVALC
*****************************************************
*
      SUBROUTINE EVALC(CVAL,IVAL,VAL,ID)
C THIS SUBROUTINE EVALUATES NUMERICAL VALUE OF CHARACTER STRING.
C BOTH INTEGER, FLOATING POINT, AND E-FORMAT NUMBERS ARE ACCEPTED.
C FOR INTEGER, RETURNS BOTH INTEGER AND REAL VALUES.
C FOR FLOATING POINT, RETURNS BOTH REAL RESULT AND ITS INTEGER PART.
C FOR E-FORMAT, RETURNS REAL RESULTS AND INTEGER EXPONENT.
C ID IS THE NUMBER OF DIGITS FOUND.  ILLEGAL CHARACTER TERMINATES
C EVALUATION.
      CHARACTER*1 V
      CHARACTER*(*) CVAL
      DATA IB /10/
      IS = ICHAR('0')
      IM = LEN(CVAL)
      NUM = 0
      NORM = 1
      IF = 0
      IE = 0
      ID = 0
      I = 0
   10 I = I + 1
      IF (I.GT.IM) GO TO 60
      V = CVAL(I:I)
   20 IV = ICHAR(V) - IS
      IF ((IV.LT.0).OR.(IV.GT.9)) GO TO 30
      NUM = IV + IB*NUM
      ID = ID + 1
      IF (IF.EQ.1) NORM=IB*NORM
      GO TO 10
   30 IF ((V.EQ.' ').OR.(V.EQ.'+')) GO TO 10
      IF (V.NE.'-') GO TO 40
      NORM = -NORM
      GO TO 10
   40 IF (IE.EQ.1) GO TO 60
      IF ((V.NE.'.').OR.(IF.EQ.1)) GO TO 50
      IF = 1
      GO TO 10
   50 IF ((V.NE.'E').AND.(V.NE.'e')) GO TO 60
      VAL = FLOAT(NUM)/FLOAT(NORM)
      NUM = 0
      NORM = 1
      IE = 1
      IF = 0
      GO TO 10
   60 IVAL = NUM/NORM
      IF (IE.EQ.0) VAL = FLOAT(NUM)/FLOAT(NORM)
      IF (IE.EQ.1) VAL = VAL*(10.**IVAL)
      RETURN
      END
*
*****************************************************
* CLEAR  -- ERASES SCREEN FOR IBM TSO
*****************************************************
*
      SUBROUTINE CLEAR
C ERASES SCREEN FOR IBM MVS/TSO
      CHARACTER*5 LBL
      DATA LBL /'CLEAR'/
C     CALL IATTCH(LBL,5,IRC,ICMDRC)
      RETURN
      END
//NAMEGEN  EXEC  PLIX,PARM='OPT(0)'
//*
//* THIS PROGRAM ADDS NAME CARDS TO OBJECT MODULES SO THAT LOAD MODULES
//* WITH SEPERATE MEMBERS FOR EACH CSECT CAN BE CREATED. WRITTEN BY CST
 /*  OBJECT MODULE NAME CARD GENERATOR  */
 NAMEGEN: PROC OPTIONS (MAIN);
 /*  THIS PROGRAM GENERATES AND INSERTS LINKAGE EDITOR NAME CARDS
     INTO A STREAM OF OBJECT MODULES GENERATED BY BATCHED COMPILATION.
     THE STATEMENT 'GET_MODULE_NAME' MAY HAVE TO BE MODIFIED TO WORK
     WITH OBJECT MODULES NOT CREATED BY THE FORTRAN COMPILERS. NOTE
     THAT, WITH THE FORTRAN COMPILERS, THE OBJECT MODULES GENERATED BY
     THE DECK OPTION AND WRITTEN TO SYSPUNCH CONTAIN SEQUENCE NUMBERS,
     WHILE THE MODULES WRITTEN TO SYSLIN DO NOT.  THEREFORE, THE DECK
     OPTION IS RECOMMENDED.  NOTE ALSO THAT NO ALIAS CARDS ARE
     GENERATED.  IF REQUIRED, ALIAS NAME CARDS MUST BE GENERATED BY
     HAND.  INPUT IS READ FROM DDNAME='INFILE', AND WRITTEN TO DDNAME=
     'OUTFILE'.
     *****
     WRITTEN BY C. THOMAS  --  08/30/74   */
 DCL
     A              CHAR (1),
     B              CHAR (3),
     C              CHAR (68),
     D              CHAR (4),
     E              FIXED BIN (31),
     F              CHAR (8),
     K              FIXED BIN (31) INIT (0);
 ON ENDFILE (INFILE)  GO TO DONE;
 LOOP:  GET SKIP FILE(INFILE) EDIT (A,B,C,D,E)  (A(1), A(3), A(68),
     A(4), F(4));
 K = K + 1;
 PUT SKIP FILE(OUTFILE) EDIT (A,B,C,D,E)  (A(1), A(3), A(68),
     A(4), P'9999');
 GET_MODULE_NAME:IF  K = 1  THEN F = SUBSTR (C,13,8);
 IF  B ^= 'END'  THEN GO TO LOOP;
 E = E + 1;
 PUT SKIP FILE(OUTFILE) EDIT ('  NAME ', F, D, E)
     (A(7), A(8), COL(73), A(4), P'9999');
 K = 0;
 GO TO LOOP;
 DONE: END NAMEGEN;
//GO.INFILE    DD  DISP=(SHR,PASS),DSN=*.COMPILE.FORT.SYSPUNCH
//GO.OUTFILE   DD  DISP=(NEW,PASS),DSN=&OBJOUT,UNIT=VIO,
//             SPACE=(TRK,(50,50),RLSE),DCB=OBJECT
//   EXEC  FORTLG,PARM.LKED='NCAL,LET'
//LKED.SYSLIN    DD  DISP=(SHR,PASS),DSN=*.NAMEGEN.GO.OUTFILE
//LKED.SYSLMOD   DD  DISP=(NEW,CATLG),
//   DSN=APP1.GRAPHICS.PCMOVIE,UNIT=___,
//   SPACE=(TRK,(50,50),RLSE),DCB=OBJECT